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1.0  Introduction 


Increased  mixing  is  sometimes  desired  in  aircraft  jet  exhaust  to  disperse  hot  gases 
quickly.  Similarly,  in  combustors  enhanced  mixing  increases  the  net  reaction  rate  and 
thereby  combustor  efficiency.  The  primary  objectives  of  this  program  were  to  investigate 
methods  for  open-loop  optimization  and  real-time  feedback  control  of  mixing  in  a  heated 
jet  and  to  demonstrate  the  resulting  techniques  in  a  laboratory  environment.  The 
techniques  developed  in  this  program  are  applicable,  for  example,  to  C-17  plume 
temperature  reduction.  The  approach  to  achieving  the  program  objective  has  involved 
systems  modeling,  control  law  design,  and  control  law  implementation  in  a  real-time 
laboratory-scale  jet. 

The  systems  modeling  involved  various  levels  of  sophistication  ranging  from  large- 
scale  direct  numerical  simulations,  which  are  described  in  Section  2.0,  to  an  inviscid 
linear  stability  approach,  which  is  described  in  Section  3.0.  The  direct  numerical  jet 
simulations,  which  accounted  for  effects  of  fluidic  actuators,  were  used  to  predict  the 
flow  dynamics  with  various  types  of  forcing.  The  mixing  effectiveness  was  determined 
using  potential  core  length  and  scalar  concentration  as  metrics.  The  direct  simulations 
were  also  used  in  combination  with  optimization  methods  to  find  the  optimal  actuation. 
In  an  attempt  to  develop  a  jet  model  which  captures  the  relevant  physics  but  has 
sufficient  simplicity  to  provide  realistic  state  estimates,  a  linear  stability  model  of  the  jet 
was  formulated.  This  approach  used  a  stability  scheme  with  a  piecewise  quadratic  base 
flow  for  obtaining  the  eigenvalues  that  provide  the  instability  development  leading  to  the 
turbulent  breakdown  of  the  jet.  It  was  found,  however,  that  the  computational  effort  for 
this  approach  made  it  unmanageable  as  a  state  model  in  a  feedback  control  setting. 

Consequently,  a  baseline  controller  was  formulated  in  the  present  work  by  using  an 
experimental  database  for  a  jet  with  pulsed  excitation.  Various  controller  design  methods 
were  evaluated  which  could  provide  the  control  logic  for  the  jet  mixing.  These  included 
a  weak  solution  of  the  linearized  flow  equations  with  a  forcing  term,  a  linearized  system 
model  determined  by  the  input-output  relationship  of  the  flow  system,  and  a  nonlinear 
inverse  system  model.  The  latter  was  used  to  construct  the  baseline  controller  and,  to 
compensate  the  inversion  error,  an  integrator  was  used  in  the  feedback  path.  The 
resulting  control  design,  which  is  described  in  Section  4.0,  was  delivered  to  GTRI  for 
laboratory  demonstration. 

The  laboratory  implementation  of  the  controller  summarized  in  Section  5.0  was 
focused  on  pulsing  air  into  the  jet  exhaust  stream  to  reduce  the  temperature  by  increasing 
mixing  of  the  exhaust  with  ambient  air.  The  miniature  engine  used  in  the  experiments  (a 
RAM  750  manufactured  by  R.A.  Microjets)  was  a  single  stage  axial  turbine  and  is 
capable  of  producing  approximately  20  lbs.  of  thrust.  Two  actuators  located  at  the  jet 
exit  pulsed  air  into  the  exhaust  stream  and  could  achieve  frequencies  up  to  1200  Hz  at 
flow  rates  up  to  10%  of  the  jet  engine  flow  rate.  This  approach  provided  the  important 
capability  of  simultaneously  achieving  both  high  frequency  and  high  flow  rate  in  a 
compact  system.  Exhaust  temperature  was  measured  at  the  immediate  jet  exit  as  well  as 
downstream  of  the  exit  with  a  movable  aspirated  thermocouple.  The  downstream 
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centerline  temperature  was  used  as  a  figure  of  merit  for  the  mixing  and  provided  the 
single  input  to  the  controller. 

In  the  GTRI  laboratory  demonstration,  all  data  acquisition  and  control  were 
performed  with  a  personal  computer  using  the  Lab  View  software  package.  The  control 
systems  consisted  of  the  jet  engine  throttle  and  the  actuator  pulsing  frequency,  while  the 
measurement  systems  consisted  of  mass  flow  to  the  actuators,  jet  exhaust  velocity,  and  jet 
exhaust  temperature.  Work  with  the  Boeing-supplied  control  algorithm,  which  was  based 
on  data  from  a  generic  axisymmetric  round  jet,  revealed  that  it  could  not  be  ported 
directly  to  a  jet  engine  which  has  both  a  centerbody  and  some  residual  swirl.  However, 
an  experimental  database  for  the  jet  exhaust  from  the  RAM  750  did  provide  the  relation 
between  centerline  temperature  and  actuation  frequency  which  allowed  a  classical 
proportional-integral  controller  to  regulate  the  plume  temperature.  This  work  is 
described  in  Section  5.0.  In  addition,  an  alternate  approach  was  adopted  in  which  a 
genetic  algorithm  was  implemented  to  provide  an  automated  means  for  optimizing 
actuator  forcing  frequency  and  amplitude.  It  was  found  that  this  technique  could  drive 
the  system  to  undesirable  states  during  optimization.  This  approach  is  also  described  in 
Section  5.0.  Conclusions  based  on  the  overall  program  effort  are  contained  in  Section 
6.0 
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2.0  Direct  Simulation  of  Turbulent  Jet  Mixing 


Numerical  simulations  have  been  used  to  study  the  dynamics  of  jets  subject  to 
different  types  of  forcing  and  to  develop  tools  for  their  optimization.  The  work  was 
divided  in  three  parts:  First,  direct  numerical  simulations  of  forced  jets  were  conducted  to 
provide  a  database  for  the  identification  of  suitable  measures/control  laws  for  jet  mixing. 
Then  numerical  simulations  were  combined  with  an  optimization  procedure  to  search  for 
the  optimal  actuation  with  respect  to  jet  mixing.  Finally,  this  procedure  was  applied  to 
the  optimization  of  jet  mixing  at  higher  Reynolds  number  using  large  eddy  simulation  of 
the  flow. 

Large-scale  direct  numerical  simulations  of  jets  forced  with  different  types  of 
actuation  were  conducted  to  perform  a  detailed  analysis  of  their  dynamics.  The 
compressible  Navier-Stokes  equations  were  solved  without  modeling  approximations. 
Fluidic  actuators  were  included  in  the  computation  by  adding  specially  designed  body 
force  terms  to  the  governing  equations.  By  taking  this  novel  approach  we  avoided  the 
computational  cost  that  would  have  been  required  to  represent  the  actuators  in  detail  but 
still  were  able  to  provide  a  high-amplitude,  low  mass  flux  forcing  that  was  similar  to  that 
used  in  corresponding  experiments.  Two  slot  jet  actuators  were  positioned  just  down 
stream  of  the  nozzle,  each  spanning  90  degrees  of  the  jet’s  circumference  and  blowing 
normal  to  shear  layer.  They  were  operated  180  degrees  out  of  phase  to  excite  flapping 
modes  in  the  jet.  When  activated,  the  actuators  mass  flux  was  3.5  percent  of  the  jet’s 
mass  flux  and  the  peak  actuator  velocity  was  60  percent  of  the  jet  velocity.  Details  of  the 
numerical  method  and  the  incorporation  of  the  fluidic  actuators  are  given  in  Appendix  A. 
Three  cases  were  simulated:  one  was  not  forced  in  order  to  provide  a  point  of  comparison 
while  the  two  others  were  forced  at  Strouhal  numbers  St  =  fD/U  =  0.2  and  St  =  0.4.  In 
the  experiments  of  Parekh  et  al.  (1996),  St  =  0.2  was  found  to  be  most  effective  at 
reducing  the  potential  core  length  and  was  picked  here  for  that  reason;  St  =  0.4  was  found 
to  be  the  most  amplified  frequency  in  the  unforced  jet.  All  three  jets  are  visualized  with 
contours  of  vorticity  magnitude  in  Figure  2-1.  Forcing  at  St  =  0.2  excited  a  large  a  large- 
scale  flapping  of  the  entire  jet  column  (Figure  2- lb)  as  observed  in  the  experiments  of 
Parekh  et  al.  When  forced  at  St  =  0.4,  the  jet  responded  more  rapidly  and  large  turbulent 
structures  appeared  closer  to  the  nozzle  (Figure  2-lc).  However,  these  instabilities 
appear  to  saturate  before  the  end  of  the  potential  core  and  the  overall  effect  on  the  jet  was 
less  pronounced  than  when  it  was  forced  at  St  =  0.2. 

Several  measures  of  mixing  effectiveness  were  investigated.  Potential  core  length 
measurements  and  scalar  concentrations  both  showed  that  forcing  at  either  Strouhal 
number  dramatically  increased  mixing.  The  streamwise  mass  flux  was  significantly 
increased  for  the  forced  cases,  with  St  =  0.2  forcing  inducing  marginally  more 
entrainment  than  St  =  0.4  forcing.  Aside  from  the  potential  core  length  difference,  the 
most  dramatic  difference  found  between  the  two  forced  cases  was  in  the  scalar 
dissipation  and  in  the  integral  of  the  squared  radial  velocity  in  the  jet.  It  was  found  that 
the  single  period  average  of  the  measures  is  very  close  to  its  longtime  average,  which  is 
important  for  the  application  of  these  control  techniques.  For  a  detailed  discussion  of  the 
jet  dynamics  and  the  various  measures  for  jet  mixing  see  Appendix  B. 
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Figure  2-1.  Vorticity  magnitude  contours  for  different  cases.  Twenty  evenly  space  contours 

between  ar0  IV j  =  0.6  and  12.0  are  shown. 

In  the  second  step,  the  simulations  were  combined  with  optimization  strategies  to 
automatically  search  for  the  optimal  actuation.  Testing  various  optimization  procedures 
it  was  found  that  evolution  strategies  are  suitable  for  the  optimization  of  this 
computationally  expensive  problem.  Starting  from  a  set  of  initial  parameter  vectors 
which  contains  the  Strouhal  numbers  and  amplitudes  of  the  actuation  x  =  (Sta,  Sti„  Aa, 
A/,),  new  vectors  were  obtained  by  random  variation.  An  initial  population  of  parameter 
vectors  is  subjected  to  repeated  mutation  and  selection  until  a  sufficiently  good  solution 
is  found.  The  advantage  of  evolution  programs  is  that  they  are  efficient,  ensuring  a  rapid 
convergence  of  the  procedure  by  evaluating  several  search  trajectories  in  parallel.  Our 
approach  for  optimizing  jet  mixing  requires  100  -  200  numerical  simulations  of  the  jet. 
In  order  to  reduce  the  CPU  time,  we  made  use  of  the  observation  that  the  jet  is  sensitive 
to  the  external  forcing  already  at  early  stages  of  the  simulation.  Since  the  flow  is 
governed  by  the  dynamics  of  the  large-scale  structures,  a  coarse  grid  was  used  for  the 
evaluation  of  the  objective  function  (Appendix  C). 
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For  the  test  case  of  an  incompressible  jet  flow  at  low  Reynolds  number  Re  =  1500,  a 
systematic  search  was  performed  for  helical  and  combined  axial  and  helical  forcing  of  the 
jet.  The  direct  numerical  simulation  used  a  second  order  finite  volume  method  to  solve 
the  incompressible  Navier-Stokes  equations  on  a  staggered  spherical  grid.  The  actuation 
was  imposed  on  the  velocity  profile  at  the  inflow.  A  spherical  coordinate  system  was 
used  which  is  able  to  follow  the  streamwise  spreading  of  the  jet  and  allows  a  well- 
balanced  resolution  of  the  flow  field  with  a  reasonable  number  of  grid  points.  An 
important  result  is  that  combined  axial  and  helical  actuation  is  much  more  efficient  with 
respect  to  jet  mixing  than  a  helical  actuation  alone.  The  combined  actuation  leads  to  a 
large  spreading  of  the  jet,  resulting  in  a  large  mixture  fraction  and  a  fast  decay  of  the 
center  line  velocity  and  scalar  concentration.  We  found  the  most  pronounced  spreading 
when  the  helical  Strouhal  was  Sth  =  0.3,  which  is  close  to  the  preferred  Strouhal  number 
of  the  jet,  and  for  an  axial  Strouhal  number  Sta «  2  '  Sth.  This  is  in  very  good  agreement 
with  experimental  results  (Parekh  et  al  (1996)).  Figure  2-2  shows  a  snapshot  of  the  scalar- 
concentration  of  the  jet  obtained  for  the  optimal  actuation  found  by  the  evolution 
strategy.  This  dual-frequency  actuation  leads  to  a  large  spreading.  The  jet  column  is 
completely  dispersed  near  the  end  of  the  domain,  indicating  good  mixing  of  the  jet  with 
the  ambient  fluid.  In  Figure  2-3  we  show  the  centerline  velocity  obtained  from  the 
flapping  and  bifurcating  jet  computations.  For  comparison  the  profile  for  a  round 
axisymmetric  jet  (standard  jet)  at  comparable  Reynolds  number  is  included.  For  both  the 
flapping  and  bifurcating  jets  the  centerline  starts  to  drop  much  earlier  than  in  the  standard 
jet.  The  decay  rate  (slope  of  the  curve  in  Figure  2-3  (right))  of  the  flapping  jet  is 
comparable  to  that  of  the  standard  jet  while  the  bifurcating  jet  decays  much  faster.  For 
the  flapping  jet  the  decay  is  approximately  linear,  for  the  bifurcating  jet  superlinear.  A 
description  of  the  best  actuation  found  for  different  types  of  forcing  and  the  resulting  jet 
flows  is  given  in  Appendix  C. 


Figure  2-2.  Snapshot  of  the  scalar  concentration  for  a  jet  actuated  with  dual-frequency 
forcing  and  St„  =  0.66,  Sth  =  0.31,  A„  =  0.025,  A*  =  0.075. 
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Figure  2-3.  The  averaged  centerline  scalar  concentration  for  the  flapping,  bifurcating  and 

standard  jets,  Re  =  1500. 

Finally,  the  optimization  procedure  was  applied  to  jets  at  higher  Reynolds  number  Re 
=  6000  using  a  large  eddy  simulation.  Including  the  dynamics  procedure  to  model  the 
small  scales  of  the  flow,  the  optimization  procedure  was  used  to  search  for  the  optimal 
actuation.  Although  the  large  eddy  simulation  of  the  jet  is  computationally  expensive,  we 
were  able  to  find  good  solutions  within  reasonable  time.  Figure  2-4  shows  a  snapshot  of 
the  jet  with  the  best  mixing  properties  obtained  for  combined  axial  and  helical  actuation. 
The  actuation  leads  to  a  large  scale  flapping  of  the  jet  column  and  a  spreading  rate,  which 
is  higher  than  that  for  flows  at  lower  Reynolds  numbers. 


Figure  2-4.  Snapshot  of  the  scalar  concentration  for  a  jet  at  Re  =  6000,  St„  =  0.79,  Sth  =  0.36, 

Aa  =  0.025  and  A„  =  0.075. 

In  conclusion,  a  detailed  study  of  forced  jets  was  performed.  Fluidic  actuators  were 
introduced  that  realistically  model  those  used  in  experiments,  and  a  tool  was  developed 
for  the  optimization  of  the  actuation.  The  measures  for  mixing  that  were  identified  by  the 
simulations  proved  to  be  suitable  for  jet  flows  subject  to  different  types  of  actuation  and 
are  promising  candidates  for  incorporation  in  a  controller. 
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3.0  Eigenmode  Modeling  of  Jet  Flows 


To  design  feedback  control  algorithms,  it  is  very  useful  to  have  a  simple  model  that 
accurately  captures  the  relevant  physics  of  turbulent  jet  flows.  Such  a  model  is 
developed  in  the  present  work  using  linear  stability  theory  to  model  the  initial 
development  of  the  instabilities  leading  to  the  turbulent  breakdown  of  a  jet.  A  key 
feature  of  the  present  work  is  a  piecewise  quadratic  approximation  of  the  mean  flow  that 
permits  rapid  solution  of  the  equations.  The  stability  equations  were  derived  by 
Alexander  Pal  in  the  present  contract  work,  and  their  formulation  is  presented  in 
Appendix  D. 

It  is  well  established  that  the  inviscid  inflectional  instability  characterized  by  Kelvin 
and  Helmholtz  dominates  the  large-scale  development  of  inflectionally-dominated  flows 
even  into  non-linear  turbulent  regimes.  An  excellent  review  discussing  many  important 
aspects  of  this  problem  is  given  by  Ho  and  Huerre  (1984).  Supporting  articles  covering 
other  aspects  of  the  problem  (including  the  nature  of  nonlinear  interaction  and  wall 
effects)  are  described  in  Cain  and  Thompson  (1986)  and  Cain,  Roos  and  Kegelmen 
(1990).  The  dominance  of  the  inviscid  inflectional  instability  in  the  present  flow 
motivates  the  use  of  the  incompressible  inviscid  stability  equations  as  the  system  model. 

A  comparison  between  the  prediction  of  linear  theory  and  the  nonlinear  evolution  of 
the  jet  shear  layers  is  given  by  Morris  et  al.  (1990)  and  Viswanathan  and  Morris  (1992) 
for  planar  and  round  jets,  respectively.  These  works  derive  systems  comprised  of 
parabolic  equations  for  the  mean  flow  development  that  depend  upon  the  local  instability 
eigenproblem  for  the  forcing.  Solution  of  this  eigensystem  is  generally  the  most 
computationally-intensive  aspect  of  generating  the  solution  to  the  evolution  model  for  the 
jet.  The  need  for  very  rapid  solutions  for  real-time  systems  modeling  motivates  the  new 
approach  to  approximate  solution  of  the  disturbance  eigenvalue  problem. 

Motivated  by  the  need  for  rapid  evaluation,  piecewise  linear  and  quadratic 
approximations  of  the  velocity  profiles  are  used  to  simplify  the  inviscid  linear 
disturbance  equations.  In  addition,  mathematical  solutions  are  somewhat  easier  to  obtain 
for  the  temporal  evolution  problem.  All  the  work  described  here  will  use  the  temporal 
analysis  combined  with  Gaster’s  relation  to  approximate  the  behavior  of  the  spatially 
evolving  problem.  The  approach  so  describing  an  inflectionally  dominated  flow  is  given 
by  Drazen  and  Reid  (1981).  A  few  examples  of  the  piecewise  linear  analysis  for  planar 
flows  will  be  presented  before  describing  a  piecewise  quadratic  approximation  that  will 
be  used  as  an  approximation  in  the  round  jet  geometry. 


3.1  The  piecewise  linear  approximation  for  planar  geometries 

The  simplest  example  of  the  piecewise  linear  approach  is  a  three-segment 
representation  of  the  plane  shear  layer.  This  approach  has  been  shown  to  provide 
consistent  and  reasonable  characterization  of  the  stability  behavior  of  this  flow,  and  it 
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compares  well  with  a  more  accurate  (but  more  time  consuming)  analysis  with  a 
hyperbolic  tangent  representation  of  the  mean  velocity  profile.  The  characterization  of 


Figure  3-1.  Phase  speed _ and  growth  rate - for  piecewise  linear  planar  shear  layer. 

the  disturbances  is  given  by  the  exponential  coefficient  for  temporal  disturbance  growth 
rate  and  phase  speed  and  is  a  function  of  the  disturbance  wavenumber.  As  shown  in 
Figure  3-1,  the  phase  speed  computed  using  such  an  approach  varies  from  U  to  3U, 
where  U  is  the  velocity  of  the  slow  stream.  The  complex  wave  speed  comes  as  complex 
conjugate  pairs,  and  the  actual  phase  velocity  is  the  mean  of  the  two  branches  plotted. 
The  imaginary  part  of  the  wave  speed  (which  in  product  with  the  wavenumber  gives  the 
exponential  grown  rate)  vanishes  at  a  wavenumber  slightly  greater  than  0.3. 

The  next  piecewise  linear  flow  considered  is  the  planar  jet.  Figure  3-2  shows  the 
Gaster  transformed  and  scaled  spatial  growth  rate  and  phase  speed  versus  wavenumber 
for  the  sinuous  mode  of  a  planar  jet  having  a  potential  core  of  18  jet  radii  (r„)  long,  a 
shear  layer  of  1  unit  width,  and  a  freestream  velocity  equal  to  1/3  of  the  jet  velocity 
1.5  Ua,  where  Ua  is  the  average  of  the  jet  and  co-flow  velocities.  Note  that  in  this  case  the 
long  wavelength  (low  wavenumber)  disturbances  have  a  phase  speed  of  the  freestream 
speed  while  shorter  disturbances  have  a  phase  speed  of  the  mean  shear  layer  speed. 
Figure  3-3  shows  the  behavior  of  the  varicose  mode.  Note  that  the  varicose  mode  has  a 
phase  speed  equal  to  that  of  the  jet  centerline  for  long  wavelength  disturbances  and  the 
same  phase  speed  as  the  sinuous  mode  for  shorter  wavelength  disturbances.  When  the 
shear  layers  are  thin  and  separated  by  a  large  region  of  potential  flow,  the  stability  of  the 
planar  jet  is  nearly  the  same  as  that  of  the  planar  shear  layer  except  at  very  long  wave 
lengths.  At  smaller  separation,  disturbances  of  the  two  shear  layers  exhibit  a  strong 
coupling.  These  behaviors  for  both  the  varicose  and  the  sinuous  modes  are  characteristic 
of  the  actual  physical  system. 
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Figure  3-2.  Phase  speed _ and  growth  rate - of  a  sinuous  mode  in  a  piecewise  linear 

planar  jet. 


Figure  3-3.  Phase  speed _ and  growth  rate - of  a  varicose  mode  in  a  piecewise  linear 

planar  jet. 


3.2  The  piecewise  quadratic  approximation  for  cylindrical  geometries 

A  piecewise  representation  of  the  mean  velocity  profile  simplifies  the  round  jet 
stability  problem  to  Bessel’s  equation.  The  solutions  are  constrained  by  requiring  finite 
levels  in  the  inner  potential  region  and  solutions  that  vanish  at  infinity  in  the  outer 
potential  region.  These  inner  and  outer  solutions  (in  terms  of  Bessel  functions)  are 
coupled  by  matching  conditions.  The  matching  is  achieved  by  a  combination  of  the  inner 
and  outer  Bessel  solutions  (a  linear  combination  of  Bessel  functions  is  a  valid  solution 
within  the  finite  thickness  shear  layer).  This  problem  was  formulated  by  Alexander  Pal 
under  the  present  contract  work  with  the  derivation  of  the  governing  equations  given  in 
Appendix  D.  In  the  present  work  the  complex-valued  quadratic  dispersion  relation  was 
solved  using  Mathematica  (Cain,  et  al.  (1998)). 

The  behavior  of  the  axisymmetric  disturbances  for  a  round  jet  with  a  potential  core 
diameter  of  18r„  and  jet  velocity  Uj  with  zero  free  stream  velocity  is  given  in  Figure  3-4. 
Note  that  the  appropriately  scaled  growth  rate  and  phase  speed  behave  in  a  manner  which 
is  qualitatively  similar  to  the  varicose  mode  of  the  analogous  planar  jet. 
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Figure  3-4.  Phase  speed  — -  and  growth  rate _ of  a  varicose  mode  in  a  piecewise  linear 

planar  jet. 

The  formulation  of  Viswanathan  and  Morris  (1992)  was  implemented  using  the 
piecewise  quadratic  stability  formulation.  An  example  of  the  predicted  evolution  of  the 
shear  layer  edges  using  only  the  n  =  0  axisymmetric  disturbance  is  shown  in  Figure  3-5. 
Figure  3-6  shows  the  locally  dominant  wavenumber  versus  downstream  distance  in  jet 
radii.  It  is  assumed  that,  when  a  disturbances  saturates  (due  to  the  thickening  of  the  shear 
layer),  the  sub-harmonic  will  become  dominant  and  evolve  until  saturation,  and  so  on. 
The  results  shown  in  Figures  3-5  and  3-6  are  for  an  initial  shear  layer  thickness  of  1%  of 
the  jet  radius,  with  the  n  =  +1  and  n  =  -1  modes  assumed  to  grow  at  the  same  rate  as  the  n 
=  0  modes. 


****•*.«, 


Figure  3-5.  Inner  and  outer  shear  layer  edges 
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Figure  3-6.  Locally  dominant  instability  wave  number. 

In  conclusion,  an  approximate  analytic  solution  to  the  appropriate  linear  stability 
problem,  and  the  role  of  this  solution  in  the  evolution  of  the  round  jet,  has  been 
investigated.  The  use  of  such  a  model  as  a  state  estimator  in  a  feedback  control 
framework  must  include  at  least  the  n  =  +  1  modes  in  addition  to  the  n  =  0  mode  given 
here.  In  the  case  of  a  thin  shear  layer,  the  n  =  ±  1  modes  behave  similarly  to  the  n  =  0 
mode  and  may  be  analyzed  with  the  n  =  0  analysis.  However,  the  solution  for  n  =  +  1  is 
required  as  the  shear  layer  thickness  becomes  significant  relative  to  the  jet  radius.  For  the 
mixing  problem  it  is  likely  that  the  solution  near  the  end  of  the  potential  core  will  be  need 
to  be  calculated.  After  further  study  it  was  concluded  that  even  the  linear  stability 
analysis  approach  was  not  manageable  as  a  state  model  in  a  feedback  control  setting. 
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4.0  Control  via  Classical  System  Identification 


The  controller  design  objective  was  to  construct  a  controller  that  regulates  the  local 
temperature  at  a  specified  downstream  location  using  a  fluidic  actuator.  For  the 
technique  to  be  implemented  in  an  aircraft,  the  controller  needs  to  operate  in  real  time. 
During  the  course  of  this  program  several  different  controller  design  methods  were 
investigated  which  are  based  on  1)  a  weak  solution  of  the  linearized  flow  equations  with 
a  forcing  term  (Abergel  and  Temam,  1990),  2)  a  linearized  system  model  determined  by 
the  input-output  relationship  of  the  flow  system  (Ikeda,  1998),  and  3)  a  nonlinear  inverse 
system  model. 

In  order  to  keep  the  computational  requirements  small  enough  to  accommodate  the 
real-time  operation,  the  first  and  second  methods  use  a  linearly  approximate  model  in 
spite  of  the  fact  that  the  actual  physical  flow  is  highly  nonlinear.  The  third  method  also 
uses  the  input-output  relationship  of  the  flow,  but  the  nonlinear  inverse  design  method  is 
used  instead  of  linearization  of  the  model.  The  resulting  controller  based  on  this 
approach  does  not  require  intense  computational  effort  in  control  update.  Therefore,  the 
nonlinear  inverse  system  model  was  adopted  to  construct  the  baseline  controller  for  the 
laboratory  implementation  and  demonstration. 

The  general  description  of  the  inverse  design  method  can  be  summarized  as  follows: 

Consider  the  system  given  by  the  first  order  ordinary  differential  equation 
x  =  f(t,x,u),  where  x,  u,  t,  respectively  denote  the  state  to  be  controlled,  admissible 
control,  and  time.  For  each  t  and  x,  /  is  assumed  to  be  a  smooth  one-to-one  function  (at 
least  of  Class  C1)  in  u.  Let  w  be  the  right  hand  side  of  the  above  differential  equation. 
That  is,  w  =f(t,x,u). 

Then  u  can  be  written  as  u  =  /'‘(t^w). 

For  the  sake  of  discussion,  it  is  assumed  that  this  inversion  is  realized  perfectly,  i.e.,  no 
approximation  and/or  rounding  errors  are  involved.  Then  by  linear  optimal  control 
theory,  one  can  construct  a  proportional  controller  gain  Kp  such  that  the  error  between  the 
command  input  and  the  feedback  states  decays  exponentially  in  time.  Hence,  the  system 
can  be  regulated.  This  can  be  seen  as  follows:  let  e(t)  be  the  error  process  given  by 

e{t)=  x{t)~  xcmd(t) 

Then  by  referring  to  the  diagram  in  Figure  4-1, 
e{t)+  K  pe{t)  =  u{t). 

By  solving  the  above  differential  equation. 


12 


e(r)  =«”*'■' e(0)+  Jje  r°u(t-o)de 


From  this  it  follows  that  the  states  can  be  driven  to  the  desired  level  by  a  positive  gain  Kp 
and  the  bounded  input. 


Figure  4-1.  Control  structure  for  the  inverse  design  method. 


In  reality,  however,  the  inversion  errors  must  be  taken  into  account.  That  is, 
u  =  f~l(t,x,w)=  w)+Au,  and  an  additional  element  is  needed  in  the  block 

diagram  for  error  correction.  This  error  correction  can  be  done  by  a  number  of  ways  (for 
example,  a  neural-network-based  on-line  learning  algorithm  or  a  simple  integrator). 

A  controller  implementation  strategy  was  then  formulated.  A  baseline  controller  was 
designed  using  an  experimental  database  for  pulsed-jet  excitation.  In  order  to  make  the 
resulting  controller  applicable  to  the  RAM  750  miniature  engine  tested  at  GTRI,  a 
normalized  temperature  T*  defined  by  the  following  relation  was  used: 

t  —  t 

T*  - —  ( where 

to~t- 

t0 :  temperature  at  the  nozzle  exit 
t„  :  ambient  temperature 
t  :  temperature  to  be  controlled 

The  baseline  controller  was  constructed  by  the  inverse  system  identification  method, 
and  an  integrator  was  used  in  the  feedback  path  to  compensate  the  inversion  error.  The 
resulting  control  architecture  is  illustrated  in  Figure  4-2. 


13 


Figure  4-2.  Control  architecture  for  active  flow  control  of  jet  temperature. 

It  should  be  noted  that  in  the  general  discussion  on  the  control  structure,  the  system 
inversion  model  was  described  by  a  first-order  differential  equation.  However,  the  actual 
baseline  controller  was  designed  by  a  zeroth  order  inverse  system  model  (i.e.,  continuous 
but  not  differential  form).  The  resulting  control  design  logic  was  coded  in  FORTRAN 
and  delivered  to  GTRI  for  laboratory  evaluation.  It  was  found  in  the  experimental 
demonstration  that  the  curve  fit  for  the  excitation  frequency  was  not  sufficiently  general 
to  describe  the  exhaust  flow  for  the  RAM  750  miniature  engine.  Consequently  it  was 
decided  to  explore  other  techniques  such  as  evolutionary  strategies  as  a  means  for 
controlling  the  temperature  in  the  jet. 
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5.0  Parameter  Optimization  via  Evolutionary  Strategies 


In  this  work  we  employ  a  new  optimization  strategy  for  the  control  of  flows  based  on 
a  class  of  algorithms  from  the  discipline  of  “soft  computing.”  Specifically,  we  consider 
the  application  of  evolutionary  strategies  for  optimization  in  both  a  simulation 
environment  and  a  real-time  physical  experiment.  The  broad  class  of  evolutionary 
strategies,  of  which  genetic  algorithms  are  a  part,  differ  from  traditional  gradient  search 
techniques  in  that  they  do  not  assume  any  knowledge  of  the  derivatives  ot  the  function 
being  optimized.  Instead,  the  optimal  (based  on  a  set  of  preferred  metrics)  value  of  a 
function  is  obtained  by  a  random  natural  selection  from  the  available  sample  set  ot 
independent  variables.  In  this  application,  the  independent  valuables  include  actuator 
number,  placement,  frequency,  phasing,  and  amplitude.  The  objective  of  the  search  in 
this  case  is  to  select  the  set  of  actuator  parameters  that  maximizes  the  jet  mixing  as 
indicated  by  plume  downstream  temperature  or  velocity.  From  iteration  to  iteration  (or 
generation  to  generation)  a  new  set  of  parameters  is  selected  based  on  a  random  selection 
process,  and  the  value  of  the  function  for  each  of  those  parameters  is  compared  against 
the  function  values  corresponding  to  a  previous  set  of  parameters.  The  parameters  that 
provide  more  optimal  function  values  survive  and  are  the  basis  for  selection  of  a  new  set 
of  parameters.  Algorithms  for  defining  the  population  of  independent  variables  that  form 
each  generation  can  be  as  simple  as  one  member  per  generation  to  numerous  members 
that  mutate  and  propagate  a  new  generation  of  members. 

The  specific  search  strategy  that  we  employ  is  a  simple  1  +  1  strategy  in  which  only 
two  function  evaluations  at  a  time  are  compared  to  determine  the  better  parameter  set. 
The  initial  evaluation  of  each  set  is  referred  to  as  the  ‘parent’  and  the  successive 
evaluation  as  the  ‘child’.  The  specific  algorithm  for  selection  of  the  next  parameter  set 
for  function  evaluation  is  illustrated  in  Figure  5-1.  This  is  akin  to  simulated  annealing 
except  that  instead  of  using  a  Boltzmann  distribution  as  the  basis  for  the  selection  of 
successive  generations,  we  employ  a  restart  process  after  a  finite  number  of  iterations  to 
assure  that  the  optimization  is  not  trapped  in  a  local  minimum. 


Algorithm 

Evaluate  f(xn) 

Evaluate  f(x*), 
x*  —  xn  +  N(o,o) 

If  f(xn)  <  f(x*)  then 
xn+1  =  xn 
else  xn+1  =  x* 

Vary  a  by  1 /5-rule 


Constraint  on  parameter  range 

Us - — - 


Figure  5-1.  Selection  strategy  for  genetic  algorithm. 
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This  approach  was  used  successfully  as  a  wrapper  around  numerical  simulations  for 
optimizing  flow  control  parameters  and  is  robust  enough  to  handle  various  types  of 
constraints  on  the  parameter  search  space.  The  notion  of  a  "wrapper”  is  simply  that  the 
genetic  algorithm  is  used  to  select  the  next  set  of  actuation  parameters  to  be  evaluated 
through  a  numerical  simulation  or  experiment.  In  one  case,  we  used  this  technique  to 
select  the  optimal  set  of  frequencies  for  maximizing  jet  mixing  with  a  DNS  simulation  as 
our  test  bed  and  in  the  other,  we  focused  on  maximizing  the  spreading  of  a  bifurcating  jet 
by  optimal  selection  of  excitation  frequency,  frequency  ratio,  and  phase  with  a  vortex 
method  simulation  (Koumoutsakos,  Freund,  and  Parekh,  1998,  included  for  reference  in 
Appendix  E). 

In  the  DNS  simulation,  the  evolutionary  algorithm  varied  the  frequency,  phase,  and 
amplitude  of  three  linearly  combined  sinusoidal  excitation  signals.  The  variation  of  these 
parameters  was  only  constrained  to  a  specified  maximum  energy  for  the  combined  signal. 
This  avoided  the  trivial  outcome  of  simply  increasing  the  amplitude  of  all  three  signals 
without  bound.  While  the  initial  guess  distributed  energy  nominally  equally  among  three 
separate  signals,  the  optimized  set  of  signals  distributed  almost  all  the  energy  into  only 
one  frequency,  a  Strouhal  number  around  0.2,  which  is  very  close  to  the  optimal  value 
found  in  previous  experiments  (see  Parekh  et  al.,  1996). 

In  the  vortex  simulations,  the  algorithm  was  able  to  “discover”  a  flow  state  not 
previously  observed.  The  bifurcating  jet  has  a  Y-shaped  structure  in  which  successive 
vortex  rings  propagate  on  alternate  trajectories  forming  two  branches.  The  angle 
separating  those  two  branches  is  indicative  of  the  degree  of  jet  spreading.  We  expected 
the  algorithm  to  select  an  optimal  set  of  parameters  that  maximized  that  spreading  angle 
but  with  no  change  to  the  general  structure  of  the  jet.  Instead,  the  algorithm  converged 
on  a  set  of  parameters  that  produced  a  bifurcating  jet  with  a  secondary  bifurcation 
resulting  in  much  greater  spreading  angles. 

These  two  cases  are  illustrative  of  the  utility  of  this  optimization  strategy  both  in 
optimizing  flow  development  without  a  priori  knowledge  about  the  flow  and  in 
discovering  new  phenomenon  that  might  easily  be  missed  even  by  systematic  searches 
through  the  parameter  space.  The  following  example  of  real-time  optimization  of  a 
small-scale  jet  engine  illustrates  the  utility  of  this  approach  in  the  presence  of  noise  and 
process  variability. 

The  laboratory  demonstration  involved  the  use  of  pulsed  fluidic  actuators  at  the 
engine  nozzle  exit  plane  for  plume  mixing  enhancement.  The  miniature  engine  used  in 
the  experiments  is  a  RAM  750  manufactured  by  R.A.  Microjets.  It  has  a  single-stage 
centrifugal  compressor  and  a  single-stage  axial  turbine.  It  is  capable  of  producing 
approximately  20  lbs.  of  thrust  and  high  subsonic  exhaust  flows  at  temperatures  over 
700°  C.  Both  the  engine  and  actuator  were  operated  under  computer  control  with  the 
genetic  algorithm  controlling  in  real  time  the  frequency  of  the  actuator  based  on 
temperature  measured  in  the  plume.  Figure  5-2  shows  the  engine  installed  on  the  test 
stand  with  the  compact  actuators  on  opposite  sides  of  the  nozzle.  The  “amplitude”  of  the 
excitation  was  set  by  controlling  the  flow  rate  through  the  actuator. 
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Three  factors  contributed  to  the  jet  characteristics  differing  significantly  from  those 
observed  in  generic  jet  flows  from  axisymmetric  laboratory  nozzles  with  top-hat  exit 
velocity  profiles.  First,  the  presence  of  the  centerbody  results  in  a  jet  exhaust  that  decays 
much  more  rapidly  than  one  would  expect  from  simply  scaling  generic  jet  data  based  on 
the  engine  nozzle  diameter.  Second,  being  a  hobby  class  engine,  the  engine  is  not 
designed  to  produce  uniform  non-swirling  flow  from  the  exit,  and  the  engine  speed 
control  tolerance  is  only  good  to  a  few  percent  due  to  it  susceptibility  to  EMI. 
Additionally,  in  these  flow  control  experiments,  one  of  the  actuators  was  turned  off  to 
simulate  a  failure  mode.  While  these  factors  are  undesirable  in  a  canonical  flow 
experiment,  they  are  ideal  for  evaluating  the  general  robustness  of  the  optimization 
technique  in  the  presence  of  realistic  complexities  and  noise  sources. 


Figure  5-2.  GTRI  laboratory  configuration  using  the  RAM  750  miniature  engine  for 
demonstration  of  turbulent  jet  mixing  control. 

A  typical  experimental  result  is  shown  in  Figure  5-3.  Starting  with  a  random  initial 
guess  for  excitation  frequency  (represented  nondimensionally  by  Strouhal  number  = 
frequency  x  nozzle  diameter  /  exit  velocity),  the  algorithm  determines  the  next  guess  by 
random  selection.  The  random  distribution  governing  the  “rolling  of  the  dice”  is  a 
uniform  distribution  of  fixed  width  centered  about  the  previous  guess.  The  function 
evaluation  consists  of  measuring  a  short-time  (up  to  2  seconds)  average  centerline 
temperature  at  a  strearnwise  location  4  diameters  downstream  of  the  jet  exit.  Whichever 
frequency  (‘parent’  or  ‘child’)  corresponds  to  a  lower  temperature  is  selected  as  the 
‘parent’  for  the  next  iteration  or  ‘generation’.  Every  ten  iterations,  the  width  of  the 
random  distribution  for  parameter  selection  is  either  increased  or  decreased  based  on  the 
1/5-rule  noted  above.  The  history  of  variations  in  terms  of  specified  excitation  frequency 
and  the  corresponding  measured  temperature  (shown  nondimensionally  as  the  difference 
between  measured  and  ambient  temperature  divided  by  the  difference  between  jet  exit 
and  ambient  temperatures)  for  a  single  experiment  is  shown  in  Figure  5-3  as  a  function  of 
iteration  number.  The  algorithm  is  able  within  a  few  iterations  to  get  very  close  to  the 
optimal  value.  As  the  optimization  converges  toward  an  optimal  value,  the  width  of  the 
random  distribution  is  automatically  scaled  down  as  is  evident  from  the  reduced 
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variations  in  the  “child”  Strouhal  numbers  at  the  latter  iterations.  The  expected 
convergence  value  of  0.10  is  different  from  the  0.12  value  to  which  the  algorithm 
converged  because  of  the  normal  variability  of  the  engine  plume  characteristics  from  run 
to  run.  This  illustrates  the  naturally  adaptive  nature  of  the  technique.  Appendix  F 
provides  additional  data  regarding  the  base  state  and  additional  optimization  experiments. 


Genetic  Algorithm  Convergence  (Ma  =  0.15,  5%  MFR,  10  Generations/cycle) 


| Parent  ST  -<—  Parent  T*  — ^ — Child  ST  — Chlld~T*~| 


Figure  5-3.  Convergence  of  excitation  frequency. 

The  final  demonstration  conducted  in  this  program  involved  implementing  the 
classical  control  teclinique  described  in  the  previous  section.  Since  the  engine 
performance  was  so  different  from  the  generic  axisymmetric  jet  used  to  develop  the  state 
model  in  the  previous  section,  a  new  empirical  model  was  derived  by  fitting  a  high-order 
polynomial  to  the  data  obtained  on  the  small-scale  engine  (see  Appendix  F).  This  model 
provided  an  analytical  approximation  of  the  relationship  between  engine  plume  centerline 
temperature  and  the  forcing  frequency  for  a  given  actuator  mass  flow  rate.  Using  a 
classical  proportional-integral  controller,  we  were  able  to  demonstrate  real-time  feedback 
control  of  plume  temperature.  Figure  5-4  shows  the  performance  of  this  controller  for 
two  values  of  relative  gain  of  the  integral  component.  During  this  experiment,  the 
commanded  temperature  was  changed  at  the  100th  iteration  of  the  control  loop.  The 
objective  of  this  particular  demonstration  was  to  show  how  real-time  feedback  control 
could  be  effected  with  even  a  very  simple  classical  control  approach.  Greater 
performance  could  certainly  be  achieved  through  optimization  or  use  of  more 
sophisticated  control  approaches.  However,  implementation  on  legacy  systems  or 
affordability  constraints  often  make  simpler  approaches  more  attractive  even  with  some 
sacrifice  in  performance. 


18 


When  real-time  feedback  control  is  required  to  track  commanded  temperatures,  the 
genetic  algorithms  are  probably  not  appropriate  since  they  could  drive  the  system  to 
undesirable  states  during  the  optimization.  A  better  role  for  them  would  be  the 
optimization  of  the  gains  and  other  parameters  of  a  more  classical  feedback  control 
technique. 


Closed  Loop  Temperature  Control  for  Two  Values  of  Integrator  Gair 
(Ma  *  0.16,  5%  MFR) 


| - Desired  Temperature  —^—Measured  Temperature  (Gain  =  1)  — * — Measured  Temperature  (Gain  =  0.25)  | 

Figure  5-4.  Real-time  feedback  control  of  jet  plume  temperature  via  a  classical 
proportional-integral  controller. 
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6.0  Conclusions 


Large-scale  direct  numerical  simulations  performed  for  jets  with  various  types  of 
actuation  proved  useful  for  optimization  of  the  actuators  with  regard  to  Strouhal  number 
and  amplitude.  The  integral  of  the  squared  velocity  and  the  integral  of  the  concentration 
proved  to  be  appropriate  objective  functions  for  evaluating  the  performance  of  a  given 
actuator.  The  simulations  revealed  that  the  spreading  of  the  jet  and  the  amount  of  mixing 
depend  strongly  on  the  actuator  Strouhal  numbers  and  that  spreading  grows  with  the 
amplitude.  Stochastic  optimization  techniques  proved  very  efficient  for  optimization  of 
the  jet  actuators,  in  particular  evolution  strategies  and  simulated  annealing.  The 
optimization  was  implemented  by  automatically  varying  the  actuation  parameters  and 
calculating  the  objective  function.  It  was  found  that  a  combined  axial  and  helical 
actuation  was  more  efficient  with  regard  to  jet  mixing  than  just  a  helical  actuation  alone. 

It  was  determined  that  a  linear  stability  model  for  a  turbulent  jet  was  not  amenable  to 
incorporation  in  a  classical  controller.  However,  the  piecewise-quadratic  representation 
of  the  base  flow  provides  a  new  computational  scheme  for  rapidly  performing  jet  stability 
analyses.  An  experimental  data  base  for  jet  mixing  based  on  a  C-17  exhaust  with 
excitation  did  not  have  the  generality  for  mixing  control  in  a  laboratory  jet  experiment 
using  a  miniature  (RAM  750)  engine.  An  experimental  database  for  the  jet  exhaust  from 
this  engine  was  used  to  provide  a  relation  between  the  centerline  temperature  and  forcing 
frequency  for  a  range  of  actuator  mass  flow  rates.  With  the  classical  proportional- 
integral  controller,  it  was  possible  to  determine  real-time  feedback  control  of  plume 
temperature.  An  evaluation  of  the  genetic  algorithm  in  the  experimental  study 
demonstrated  that  this  approach  could  drive  the  system  to  undesirable  states  during  the 
optimization.  However,  it  would  be  possible  to  use  genetic  algorithms  to  optimize  the 
gains  and  other  parameters  associated  with  the  classical  feedback  control  technique. 
Such  an  effort  would  be  a  topic  for  follow-on  research  in  this  area. 
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Appendix  A 

Mixing  enhancement  in  jet  exhaust  using  fluidic  actuators:  direct 

numerical  simulations 


NOMENCLATURE 

a  Speed  of  sound 
D  Jet  diameter 
e  Fluid  total  internal  energy 
r  Jet  radial  coordinate 
r0  Jet  radius 

Sto  Strouhal  number  =  fD/Uj 
Ta  Actuator  fluid  temperature 
Tj  Jet  exhaust  fluid  temperature 
Ambient  fluid  temperature 
Ua  Peak  actuator  fluid  velocity 
Uj  Jet  exit  centerline  velocity 
ux  Axial  velocity 
vr  Radial  velocity 
v$  Azimuthal  velocity 

x  Jet  axial  coordinate  measured  from  nozzle 
9  Jet  azimuthal  coordinate 
p  Fluid  density 
£  Scalar  concentration 

INTRODUCTION 

This  work  is  part  of  a  larger  project  which  is  exploring  means  to  control  free  shear  flows 
to  enhance  mixing  or  reduce  acoustic  emission.  Presently,  these  are  separate  goals  because 
increased  mixing  will  often  increase  noise.  Increase  mixing  is  sometime  desired  in  aircraft  jet 
exhaust  to  disperse  hot  gases  quickly.  Similarly,  in  combustors  enhanced  mixing  increases 
the  net  reaction  rate  and  thereby  combustor  efficiency.  Suppression  of  acoustic  emission  has 
obvious  application  in  the  civilian  aircraft  industry  where  noise  reduction  is  a  key  design  goal 
for  both  present  day  subsonic  and  proposed  supersonic  aircraft.  Presently,  we  distinguish 
mixing  augmentation  and  noise  reduction  as  separate  goals  because  enhanced  mixing  will 
often  increase  noise. 

Previous  studies  of  jet  forcing  have  been  conducted  in  laboratory  experiments.  Parekh  et 
al.  (1996)  excited  the  jet  shear  layer  adjacent  to  the  jet  nozzle  exit  with  both  high  frequency 
piezoelectric  actuators  and  fluidic  actuators.  We  will  focus  on  the  fluidic  actuators  which 
act  on  the  jet  by  blowing  fluid  toward  or  pulling  fluid  away  from  the  shear  layer.  Fluidic 
actuators,  particularly  those  that  are  purely  blowing,  show  promise  for  application  in  real 
jet  engines  because  of  the  availability  of  high  pressure  in  the  compressor.  Naturally  the 
‘cost’  of  control  will  depend  upon  how  much  mass  flow  and  at  how  much  pressure  is  needed 
to  effect  the  desired  response  in  the  jet.  In  small  scale  experiments,  Parekh  et  al.  (1996) 
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found  that  remarkably  little  mass  flow  (1%)  could  dramatically  increase  the  mixing  and 
change  the  overall  behavior  of  the  jet  downstream.  In  one  configuration  they  placed  thin 
slot  jet  actuators  blowing  normal  to  the  shear  layer  on  opposite  sides  of  the  jet  (see  Fig.  .1). 
When  forced  at  a  Strouhal  number,  Sto  =  0.2  and  with  mean  injection  velocity  of  0.27 Uj, 
the  length  of  the  potential  core  was  reduced  from  4  to  2  jet  diameters.  The  jet,  which  was 
axisymmetric  in  the  mean  without  forcing,  was  seen  to  flap  violently  in  the  plane  of  the 
actuators.  This  was  reflected  in  the  mean  velocity  profiles  downstream  which  became  very 
broad  in  the  plane  of  the  actuators. 

It  is  desirable  from  a  control  theory  perspective  to  develop  a  simple  model  for  the  jet  and 
its  response  to  forcing.  Linear  analysis  may  offer  such  a  model,  but  its  applicability  is 
uncertain  in  this  case  because  the  actuator  forcing  amplitude  is  quite  high.  For  example, 
the  forcing  used  by  Parekh  et  al.  (1996)  has  a  root-mean-squared  value  of  the  same  order  as 
the  turbulent  kinetic  energy  in  the  absence  of  forcing.  The  suggestion  that  linear  methods 
may  be  applicable  is  due  to  the  remarkable  success  they  have  seen  in  unforced  free  shear 
flows  ( e.g .  Morris  et  al.  1990)  in  the  presence  of  high  amplitude  disturbances.  In  the 
language  of  control  theory,  direct  simulation  of  the  complete  system  offers  an  excellent 
means  of  evaluating  estimates  for  its  state  and  response  to  control.  Direct  simulations  also 
offer  a  detailed  description  of  the  dynamics  which  is  unavailable  by  other  means. 

The  present  study  explores  the  possibility  of  direct  numerical  simulation  of  forced  jets 
and  develops  appropriate  techniques  for  their  computation.  There  are  several  key  issues 
that  need  to  be  addressed.  Boundary  conditions,  particularly  those  capable  of  providing 
very  thin  shear  layers  near  the  nozzle  exit,  are  important.  If  the  shear  layers  are  too 
thick,  the  receptivity  of  the  layer  will  change  dramatically  and  so  will  the  response  of 
the  jet  to  the  forcing.  Incorporation  of  the  actuators  themselves  is  also  important  and 
must  be  accomplished  in  a  robust  fashion  that  does  not  compromise  their  effects.  From  a 
practical  standpoint,  it  is  also  necessary  that  the  inclusion  of  the  actuators  not  require  a 
large  number  of  mesh  points.  Solutions  to  these  problems,  along  with  the  particulars  of  the 
flow  are  discussed  in  the  following  section.  This  is  followed  by  a  discussion  of  some  results, 
including  a  comparison  with  experiments. 


PRELIMINARIES 

Flow  description 

The  subject  of  the  present  study  is  a  Mach  0.8  round  jet  exhausting  into  quiescent  ambient 
fluid.  The  jet  temperature  ratio  was  Tj /T^  =  2.  Fluidic  actuators  in  the  form  of  slot  jets 
were  located  on  opposite  sides  of  the  initial  jet  shear  layer  near  the  nozzle  exit  (see  Fig.  .1). 
This  geometry  is  similar  to  that  of  Parekh  et  al.  (1996).  The  temperature  of  the  actuator 
fluid  was  Ta/Too  —  1.5.  The  actuators  were  forced  out  of  phase  with  each  other  at  Stjj  =  0.2 
with  exit  velocities  sinusoidally  varying  between  0  and  0.35 Uj.  The  jet  and  actuator  fluids 
were  marked  separately  with  passive  scalars  (£i  and  £2  respectively)  which  diffused  with 
Schmidt  number  Sc  =  1.  The  Prandtl  number  was  Pr  =  0.7  and  the  Reynolds  number 
was  limited  to  Rej  =  pjUjDj/pj  =  800.  This  is  significantly  lower  than  in  the  experiments 
of  Parekh  et  al.  and  orders  of  magnitude  smaller  than  in  actual  jet  engines,  but  it  is  well 
known  that  many  global  characteristics  of  free  shear  flows  are  independent  of  Reynolds 
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Figure  .1:  Geometry  schematic  showing  the  jet  nozzle  (a)  and  the  actuators  (b). 

number.  In  jets  at  similar  Mach  number,  it  has  been  found  that  decreasing  the  Reynolds 
number  tends  to  increase  the  potential  core  length  but  otherwise  the  mean  velocity  profiles 
are  very  similar  at  different  Reynolds  numbers  (Stromberg  et  al.  1980).  At  low  Reynolds 
numbers,  the  jet  shear  layer  takes  longer  to  transition  to  turbulence  (this  causes  the  longer 
potential  cores),  but  this  does  not  seem  to  effect  the  present  study  where  the  forcing  excites 
a  single  mode  which  then  dominates  the  flow  above  any  turbulent  fluctuations. 

Numerical  Methods 

The  fully  compressible  Navier-Stokes  equations  and  two  passive  scalar  transport  equations 
were  solved  in  cylindrical  coordinates  without  modeling  assumptions.  Here  we  give  only 
a  summary  of  the  method,  details  are  given  in  Freund  et  al.  (1997)  where  the  same  code 
was  used  to  study  turbulence  in  compressible  mixing  layers.  Sixth  order  compact  finite 
differences  (Lele  1992)  were  used  in  the  streamwise  ( x )  and  radial  (r)  directions  and  a 
Fourier  spectral  method  was  used  in  the  azimuthal  direction  (0).  Time  advancement  was 
accomplished  with  an  explicit  fourth  order  Runge  Kutta  algorithm.  At  the  r  =  0  coordinate 
singularity,  the  equations  were  solved  in  Cartesian  coordinates.  To  maintain  a  reasonable 
numerical  time  step  given  the  restriction  imposed  by  the  Runge  Kutta  algorithm,  higher 
Fourier  modes  in  6  were  not  used  near  r  =  0.  They  were  omitted  systematically  so  that 
the  effective  azimuthal  resolution  remained  nearly  constant  with  radial  location.  The  com¬ 
putational  mesh  was  320  x  180  x  96  points  in  the  axial,  radial,  and  azimuthal  directions 
respectively  and  points  were  compressed  in  the  radial  direction  around  r  =  r0. 

The  boundary  conditions,  the  jet  nozzle,  and  the  actuators  were  all  included  in  the  calcula¬ 
tion  by  modifying  the  equations  in  localized  regions  of  the  computational  domain.  For  the 
boundary  conditions,  terms  were  added  to  the  equations  in  a  zone  surrounding  the  physical 
portion  of  the  domain.  The  purpose  of  these  terms  was  to  reduce  the  effect  of  truncating 
the  computational  box  at  a  finite  distance  and  to  form  an  appropriately  thin  jet  shear  layer 
as  an  inflow  condition.  Taking  N(q)  =  0  to  be  the  governing  equations  in  conservative 
form,  where  q  =  [p,  pu,  pvr,  pvg,  e,p^i,p6]T  is  a  vector  of  the  conservative  variables,  the 
modified  equations  are 

N(q)  =  -,r(q-qg)-Cf,|2  (.1) 

where  qg  is  the  quiescent  target  state  to  which  the  solution  is  driven.  These  terms  are 
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Figure  .2:  Computational  domain  schematic  showing  the  “physical  domain”  and  the  bound¬ 
ary  zones:  (a)  inflow  zone,  (b)  far-field  zone,  and  (c)  outflow  zone.  The  small  boxes  indicate 
the  approximate  location  of  the  actuators. 


non-zero  only  in  the  boundary  zones  shown  schematically  in  Fig.  .2.  The  physical  portion 
of  the  domain  was  25.5r0  long  and  extended  to  6 r0  in  the  radial  direction.  Zone  (a)  shown 
in  Fig.  .2  was  2 r0  wide,  zone  (b)  was  1.3r0  wide,  and  zone  (c)  was  2.5r0  wide. 

The  flow  was  driven  toward  a  thin  jet  shear  layer  in  zone  (a)  to  mimic  the  effect  of  a  nozzle. 
The  value  of  a  in  this  zone  was  2a00/r0  and  the  target  axial  velocity  flow  was 


ux 


=  f(l-tanh[i(r-i) 


(•2) 


with  (3  =  0.06.  The  density  and  scalar  profiles  were  similarly  specified.  Taking  x  =  uXg /Uj, 
the  density  profile  was  pg  =  (1  —  x)Poo  +  XPj  an<i  the  scalar  profile  was  £ig  =  pgX-  These 
profiles  correspond  to  momentum  thickness  of  approximately  5m  —  0.02ro.  The  complete 
target  condition  was  then  qg  =  [pg,  pguIg,  0, 0, eg,pg£ig, 0]r,  were  eg  is  the  total  energy.  To 
calculate  eg,  we  assume  that  the  static  pressure  at  the  nozzle  exit  is  equal  to  the  ambient 
pressure.  Away  from  the  ‘nozzle’,  the  inflow  non- reflecting  boundary  condition  of  Freund 
(1997)  was  applied. 

In  the  outer  radial  boundary  zone  (b),  a  was 


(r_  -  6r0  \  3  Qqo 
\  1.3r0  )  r 


and  qg  was  set  to  the  ambient  conditions.  The  outflow  boundary  condition  (c)  was  exactly  as 
specified  in  Freund  (1997).  In  (c),  qg  was  estimated  by  assuming  self-similar  incompressible 
flow  (the  centerline  Mach  number  is  significantly  reduced  as  the  right  side  of  the  domain 
is  approached)  and  extrapolating  from  inside  the  domain.  The  flow  well  upstream  of  the 
boundary  was  found  to  be  insensitive  to  the  particular  outflow  boundary  conditions.  It  will 
be  seen  that  the  flow  changes  dramatically  when  forced,  and  work  is  currently  underway 
assessing  how  best  to  adjust  the  boundary  conditions  for  this. 

To  provide  a  seed  for  the  turbulence,  random  velocity  fluctuations  were  added  to  the  flow. 
This  was  done  with  a  cr-term  as  in  (.1)  acting  at  the  first  mesh  point  inside  the  inflow 
boundary.  We  found  that  a  =  2.8  and  qg  =  [p, p^md) P^rnd, P'/'rnd,  e, p£i,p&]T,  where  </>rnd 
are  random  numbers  <prn a  €  (-1, 1),  provided  a  sufficient  seed  for  the  turbulence.  Each 
mesh  point  and  velocity  component  was  assigned  a  new  random  number  every  time  step. 
The  actuators  extended  from  x\  =  0.26ro  to  X2  =  0.92 r0  in  the  axial  direction  and  from 
r\  =  1.24r0  to  r<2  =  3.51r0  in  the  radial  direction.  In  this  region,  a  =  2aO0/r0  and  qg  = 
[•75poo,0, -.75poo^a)0,  ea,0,  .75poo]T  with  ea  calculated  from  .75poo,  va,  and  the  ambient 
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Figure  .3:  The  x-r  shape  of  the  actuation  velocity. 

pressure.  Pure  $2  =  1  marked  the  actuator  fluid.  The  instantaneous  actuator  target  velocity 
was 

va  =  Uafx(x)fr{r)fe{0)ft(t).  (-4) 

Its  amplitude  was  set  at  Ua  =  A4Uj,  but  the  actual  peak  vr  at  the  actuator  typically  only 
reached  .35 Uj.  The  spatial  and  temporal  dependencies  of  va  were 

fx(x)  =  e-M*-o.5(xi+a:2))2  (.5) 

fr(r)  =e-a,(r-n)2  (.6) 

m  =  exp[-C(0)8e"<w2+2-8]  +  exp[-(C(0)  -  4)8e~KW-4)2+2-8]  (.7) 

ft(t)  =  0.5[1  +  sin(2ir5tDl7ji/D)sgn(sin0)],  (.8) 

where  ax  -  20.4,  ar  =  1.34,  and  ((9)  =  S9/2n  -  2.  The  x-r  shape  of  va  is  plotted  in  Fig. 
.3  and  the  azimuthal  dependence  is  plotted  in  Fig.  .4.  The  ‘sgn’  term  in  (.8)  makes  the 
actuators  act  out  of  phase. 

RESULTS 

The  velocity  field  in  a  x-r  plane  region  around  an  actuator  near  peak  blowing  is  shown  in 
Fig.  .5.  To  make  viewing  of  the  vectors  easier,  only  every  other  mesh  point  in  the  axial 
direction  and  every  fourth  mesh  point  in  the  radial  direction  are  shown.  There  are  only  9 
mesh  points  across  the  actuator  in  the  streamwise  direction  and  those  in  the  middle  of  the 
actuator,  where  the  Gaussian  axial  distribution  peaks,  carry  most  of  the  mass  flux.  The 
vectors  show  that  fluid  is  blown  out  of  the  actuator  and  into  the  jet  shear  layer  with  a 
pattern  that  resembles  an  asymmetric  stagnation  point  on  a  wall.  Also  like  a  stagnation- 
point  flow,  there  is  a  region  of  elevated  pressure  near  the  stagnation  point  just  below  the 
actuator.  This  causes  the  main  jet  flow  to  deflect  slightly  downward.  The  bulk  of  the 
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Figure  .4:  The  9  shape  of  the  actuation  velocity. 


actuator  fluid  is  turned  downstream  by  the  main  jet  flow,  but  a  small  amount  also  turns 
upstream.  This  will  also  be  shown  below  where  a  passive  scalar  is  used  to  track  the  actuator 
fluid. 

The  downstream  development  of  the  unforced  and  forced  shear  layers  are  shown  in  Figs.  .6 
and  .7  respectively.  In  the  unforced  case,  the  shear  layer  spreads  slowly  downstream  and 
it  eventually  pinches  off  the  potential  core  at  around  18  jet  radii  from  the  ‘nozzle’  exit. 
Instabilities  in  the  mixing  layers  grow  slowly  and,  though  structures  are  apparent,  no  mode 
appears  dominant.  In  the  forced  case,  the  single  mode  forced  by  the  actuators  dominates 
the  flow.  After  the  potential  core  closes,  the  jet  spreads  very  quickly  in  the  plane  containing 
the  actuators,  while  in  the  plane  perpendicular  to  it  the  jet  spreads  slowly.  The  rapidly 
spreading  portion  of  the  jet  reaches  the  boundaries  of  the  computational  domain  in  only  a 
short  distance. 

It  is  of  interest  to  track  the  jet  and  actuator  fluids  with  passive  scalars.  Figure  .8  shows 
concentration  contours  of  both  jet  and  actuator  fluids.  The  actuator  fluid  is  pulled  down¬ 
stream  by  the  jet  flow  and  quickly  mixes  into  the  jet  shear  layer.  Again  there  is  evidence 
that  a  small  amount  is  also  sent  upstream  toward  the  ‘nozzle’,  but  the  nearby  pretense  of 
the  inflow  boundary  zone  prevents  on  from  reaching  any  conclusions  regarding  how  far  up¬ 
stream  it  does  travel.  By  the  point  where  the  potential  core  closes,  the  actuator  fluid  is  very 
diffuse  and  its  mixture  fraction  is  below  £2  =  0.05  at  all  points  downstream.  Contours  of 
jet  fluid  concentration  show  evidence  of  large  scale  structures  near  the  end  of  the  potential 
core  which  send  fluid  in  mushroom  like  ejections  away  from  the  centerline  of  the  domain. 
Jet  fluid  clearly  accounts  for  the  bulk  of  the  flow  downstream  and  the  ratio  of  jet  fluid  to 
actuator  fluid  in  the  whole  domain  is  30:1.  Though  their  mass  flux  is  small,  these  actuators 
have  significantly  higher  mass  flux  relative  to  the  jet  than  do  the  experiments  of  Parekh  et 
al.  (1996).  It  is  desirable  to  make  the  actuators  smaller  and  the  jet  shear  layers  thinner  to  be 
in  better  accord  with  those  experiments,  but  to  do  this  would  require  more  mesh  points.  To 
reduce  the  size  of  the  actuators,  it  may  also  be  necessary  to  increase  the  Reynolds  number 
so  that  viscosity  does  not  overwhelm  their  action.  Increasing  the  Reynolds  number  would 
also  increase  the  resolution  requirements  for  the  entire  turbulent  portion  of  the  flow  and 
this  may  quickly  become  prohibitively  expensive  with  direct  numerical  simulation.  So,  the 
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Figure  .5:  Actuator  velocity  field  at  0  =  tt/2.  The  boxed  region  shows  region  of  non-zero 
forcing.  Contours  show  pressure  at  (p  —  Poo)  /  PjUj  =  0.006,0.043  and  0.079. 


Figure  .6:  Vorticity  magnitude  contours  for  the  unforced  case.  There  are  40  evenly  spaced 
contours  with  maximum  cor0/uj  =  8.43. 
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Figure  .7:  Vorticity  magnitude  contours  for  the  forced  case.  Rectangles  indication  the 
actuator  locations.  There  are  40  evenly  spaced  contours  with  maximum  u)r0/u3  =  8.43. 


effect  of  reducing  the  actuator  mass  flow  while  maintaining  the  peak  actuator  velocity  is 
not  known.  However,  it  was  straightforward  to  increase  the  actuator  mass  flux.  Increasing 
it  by  roughly  a  factor  of  3  only  further  reduced  the  potential  core  length  slightly  (~  15%) . 
Though  these  simulations  have  only  run  for  a  short  time,  the  statistical  sample  was  sufficient 
to  estimate  the  mean  axial  velocity.  Figure  .9(a)  shows  the  centerline  velocity  for  the  forced 
and  unforced  cases.  The  curves  are  not  smooth  because  of  the  small  statistical  sample 
available.  In  the  forced  case,  the  centerline  profile  drops  abruptly  a  short  distance  from  the 
nozzle.  If  we  define  the  end  of  the  potential  core  as  the  point  where  u (r  —  0)  starts  to  drop, 
then  we  see  that  its  length  is  halved  by  the  forcing.  These  results  are  very  similar  to  those 
of  Parekh  et  al.  (1996)  who  also  observed  a  factor  of  2  drop  in  the  potential  core  length. 
Velocity  profiles  at  19  radii  downstream  from  the  nozzle  are  shown  in  Fig.  .9(b).  The  peak 
velocity  drops  dramatically  between  the  forced  and  unforced  cases.  Surprisingly,  the  mean 
flow  becomes  nearly  uniform  across  the  entire  computational  domain  in  the  9  =  7t/2  plane. 
Similarly  flat  profiles  were  observed  by  Parekh  et  al.  The  profile  in  the  9  —  0  plane  has  the 
same  Gaussian  shape  as  the  unforced  case.  Reversed  flow  is  observed  for  the  forced  case  in 
the  9  =  0  plane.  This  is  believed  to  be  either  a  transient  associated  with  the  startup  of  the 
actuators  or  an  effect  of  the  boundary  conditions.  Work  is  underway  to  resolve  this. 

The  goal  of  the  actuation  was  to  increase  downstream  mixing.  To  quantify  this  further,  we 
will  estimate  the  total  amount  of  heat  radiation  from  the  jet.  If  we  assume  that  every  point 
acts  as  a  black  body  radiator,  then  the  total  radiated  energy  is  proportional  to 

/  (T4  -  T^)dV  (.9) 

Jv 

where  take  V  is  the  physical  portion  of  the  computational  domain.  Using  this  very  simple 
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Figure  .8:  Scalar  concentration.  The  dark  contours  show  actuator  fluid  with  levels  0.05  to 
0.25  and  spacing  0.05  and  level  1.0  is  pure  actuator  fluid.  The  light  contours  show  jet  fluid 
with  levels  0.05  to  0.95  with  0.05  spacing  and  1.0  is  pure  jet  fluid.  The  ‘nozzle’  location  is 
indicated  by  the  heavy  dashed  lines. 


(a) 


0  2  4  6 


x/r0 

Figure  .9:  (a)  Jet  centerline  velocity  for  the  unforced - and  forced - cases,  (b)  Mean 

axial  velocity  at  x  =  19r0  for  the  unforced  jet . ;  the  forced  jet  in  the  0  =  0  plane - ; 

and  the  forced  jet  in  the  6  =  tt/2  plane - . 
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estimate,  the  total  radiated  energy  drops  by  a  factor  of  2  due  to  the  actuation. 

CONCLUSIONS 

Direct  numerical  simulations  of  forced  jets  have  been  performed  which  match  qualitatively 
with  the  experiments  of  Parekh  et  al.  (1996).  New  techniques  for  generating  thin  shear  layer 
jet  inflow  boundary  conditions  and  for  including  small  actuators  into  the  computations  were 
successfully  developed.  It  was  found  that  the  actuators  need  only  slightly  perturb  the  thin 
jet  shear  layer  near  the  nozzle  to  cause  a  dramatic  change  in  the  overall  flow  development 
downstream.  Mixing,  as  measured  by  an  estimate  of  the  total  energy  radiation  from  hot 
jet  gases,  was  suppressed  by  a  factor  of  2  due  to  the  actuation.  The  potential  core  was 
similarly  shortened  by  a  factor  of  2. 

Work  continues  to  study  aspects  of  the  forcing  which  are  difficult  to  analyze  experimentally 
and  a  project  is  presently  underway  to  calculate  the  optimum  temporal  profile  for  forcing.  It 
will  also  be  useful  to  test  how  well  linear  modal  analyses  may  predict  the  observed  reaction 
of  the  jet  to  forcing.  Acoustics  are  difficult  to  study  experimentally  because  it  is  hard  to 
simultaneously  measure  the  acoustic  sources  and  the  sound  field.  Simulations  may  offer 
insight  into  means  of  reducing  acoustic  emission  by  forcing  near  the  nozzle. 
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Jet  mixing  enhancement  by  high  amplitude  fluidic  actuation 


Summary 

Recent  experiments  have  shown  that  properly  designed  high  amplitude,  low  mass  flux  pulsed  slot 
jets  blowing  normal  to  a  jet’s  shear  layer  near  the  nozzle  can  significantly  alter  the  jet’s  development 
(Parekh  et  al.  AIAA  Paper  96-0308).  In  contrast  to  commonly  used  low  amplitude  forcing,  this 
strong  excitation  appears  to  overwhelm  the  turbulence,  having  nearly  the  same  effect  at  high  and 
low  Reynolds  numbers.  It  can  therefore  be  studied  in  detail  by  direct  numerical  simulation.  In 
this  study,  direct  numerical  simulations  of  Mach  0.9,  Reynolds  number  3600  jets  exhausting  into 
quiescent  fluid  are  conducted.  Physically  realistic  slot  jet  actuators  are  included  in  the  simulation 
by  adding  localized  body- “force”  terms  to  the  governing  equations.  Three  cases  are  considered  in 
detail:  a  baseline  unforced  case  and  two  cases  that  are  forced  with  flapping  modes  at  Strouhal 
numbers  0.2  and  0.4.  (St  —  0.4  was  found  to  be  the  most  amplified  frequency  in  the  unforced 
case.)  Forcing  at  either  frequency  causes  the  jet  to  expand  rapidly  in  the  plane  parallel  with  the 
actuators  and  to  contract  in  the  plane  perpendicular  to  the  actuators,  as  observed  experimentally. 
It  is  found  that  the  jet  responds  closer  to  the  nozzle  when  forced  at  St  =  0.4,  but  forcing  at 
St  =  0.2  is  more  effective  at  spreading  the  jet  further  downstream.  Several  different  measures  of 
mixing  (scalar  dissipation,  volume  integrals  of  jet  fluid  mixture  fraction,  and  point  measurements 
of  mixture  fraction)  are  considered,  and  it  is  shown  that  by  most,  but  not  all,  measures  forcing  at 
St  =  0.2  is  the  more  effective  of  the  two  at  mixing. 


Nomenclature 

a  =  speed  of  sound 

D  =  jet  diameter 

V  =  planar  integral  of  |V£| 

e  =  total  energy 

J  —  volume  integral  of  £n 

M.  =  planar  integral  of 

r  =  radial  coordinate 

ra  =  jet  nozzle  radius 

St  =  forcing  Strouhal  number  =  fD/Uj 

T  —  temperature 

Ua  =  peak  actuator  fluid  velocity 

Uj  —  jet  exit  velocity 

vx  =  axial  velocity 

vr  =  radial  velocity 

v$  =  azimuthal  velocity 

x  =  axial  coordinate 

6  =  azimuthal  coordinate 

p  =  fluid  density 

£  =  jet  fluid  mixture  fraction 

uj  =  vorticity  magnitude 
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Introduction 


There  are  several  technological  applications  where  enhanced  jet  mixing  can  lead  to  improved  effi¬ 
ciency,  reliability,  or  safety.  For  example,  enhanced  jet  mixing  can  reduce  temperatures  on  in-plume 
aerodynamic  surfaces,  such  as  the  blown  flap  on  a  C-17  aircraft,  and  thus  provide  greater  flexibility 
in  the  choice  of  materials  for  their  construction.  Similarly,  the  mixing  efficiency  of  fuel  jets  in  com¬ 
bustors  is  an  important  factor  in  their  overall  performance,  with  size  and  weight  reductions  possible 
if  mixing  is  improved.  In  the  present  work  we  focus  on  free  jets,  with  the  principle  objective  being 
plume  temperature  reduction. 

In  general,  attempts  to  control  jets  can  be  divided  into  two  categories:  active  and  passive. 
Examples  of  passive  control  are  tabs  located  at  the  nozzle  exit,1,2  crown  shaped  nozzles,3  and 
various  other  tailorings  of  the  nozzle  exit.4-6  This  list  is  far  from  exhaustive.  Passive  control  is 
attractive  because  in  many  cases  it  entails  only  simple  design  modifications,  a  change  in  nozzle 
geometry  for  example.  Also  its  simplicity  typically  makes  the  resulting  hardware  less  subject  to 
failure.  However,  active  control,  where  nozzle  conditions  are  continuously  updated,  has  greater 
flexibility  and  therefore  greater  potential  to  modify  the  jet  flow.  In  this  study  we  analyze  a  recently 
proposed  method  for  active  control  of  high  Reynolds  number  jets.7 

In  the  past,  active  excitation  has  been  used  extensively  to  understand  the  dynamics  of  free 
shear  flows,  particularly  the  dynamics  of  the  largest  turbulent  flow  structures.  Studies  up  to  the 
mid  1980’s  are  summarize  by  Ho  &  Huerre8  and  relatively  more  recent  efforts  are  discussed  by  Koch 
et  al .9  and  Ho  et  al.10  Here,  the  excitation  used  was  typically  low  amplitude,  often  serving  only  to 
seed  instability  waves  in  the  flow  in  order  to  phase  correlate  coherent  structures.  Unfortunately, 
it  seems  that  to  modify  the  flow  significantly  at  high  Reynolds  numbers,  low  amplitude  forcing  is 
ineffective  because  the  applied  perturbations  are  overwhelmed  by  the  turbulence.  Thus,  control  by 
low  amplitude  excitation  is  not  practical  in  many  flows  of  engineering  interest. 

Recently,  a  scheme  has  been  developed  to  control  high  Reynolds  number  jets,  such  as  the 
exhaust  flow  from  jet  engines,  by  forcing  with  actuation  velocities  greater  than  the  local  turbulence 
intensity.  This  approach  was  tested  by  Parekh  et  al.7  who  designed  slot  jets  that  blew  normal  to 
the  jet  shear  layer  adjacent  to  the  nozzle  (as  in  figure  0.1).  When  pulsed  180°  out  of  phase  from  one 
another  with  peak  blowing  at  approximately  one  third  the  jet  velocity  (and  approximately  2  percent 
of  the  jet  mass  flux),  they  excited  large-scale  oscillations  in  the  jet  that  reduced  the  potential  core 
length  by  over  a  factor  of  two.  Recent  test  results  have  shown  mixing  enhancement  using  this 
technique  on  a  full-scale  engine.11  Similar  results  have  been  obtained  using  zero  net  mass  flux 
“synthetic  jets”  .7’ 12, 13  It  appears  that  these  active  control  approaches  have  been  more  successful 
at  increasing  mixing  at  high  Reynolds  numbers  than  any  attempts  by  passive  control  approaches. 
Obviously  the  direct  numerical  simulations  used  in  this  effort  are  incapable  of  addressing  high 
Reynolds  number  flows,  but  we  shall  see  that  the  behavior  of  the  simulated  forced  jets  is  similar  to 
that  observed  at  high  Reynolds  numbers,  and  new  insights  are  provided. 

Thus  far,  the  forcing  parameters  in  these  experiments  have  been  picked  beforehand,  as  in 
an  “open-loop”  control  strategy.  “Closed-loop”  control,  where  the  jet  flow  would  be  continuously 
monitored  and  its  state  used  to  update  the  control  parameters,  may  offer  improved  performance.  To 
implement  this  approach  a  practical  measure  of  the  performance  (an  objective  or  “cost”  function) 
is  needed.  An  objective  of  the  present  effort  is  to  study  metrics  for  mixing  and  provide  a  database 
for  direct  comparison  of  different  metrics.  The  choice  of  a  metric  for  a  particular  application  will 
depend  not  only  upon  its  relevance  for  the  given  mixing  objective  but  also  on  practical  aspects  of 
its  implementation.  Here  we  concentrate  only  upon  the  metrics  themselves. 
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Figure  0.1:  Geometry  schematic  showing  the  jet  nozzle  (a)  and  the  actuators  (b). 


Flow  parameters 

The  focus  of  this  paper  is  a  round  jet  at  Mach  0.9.  The  jet  Reynolds  number,  based  upon  centerline 
flow  conditions  at  the  nozzle  exit,  is  3600,  and  the  stagnation  temperature  of  the  jet  is  constant. 
These  parameters  match  those  studied  experimentally  by  Stromberg  et  al.1A  and  a  baseline,  unforced 
case  has  been  validated  against  their  data.15  Direct  comparison  shows  that  mean  Mach  number 
profiles  (and  the  overall  sound  pressure  levels  on  an  arc  at  60  radii  from  the  nozzle)  agree  to  within  5 
percent.  Unfortunately,  the  precise  nozzle  conditions  were  not  measured,  but  spectra  show  that  the 
initial  jet  shear  layers  were  laminar,  as  expected  at  this  Reynolds  number.  To  model  appropriate 
nozzle  conditions  a  rounded  ‘top-hat’  velocity  profile  was  specified  in  the  present  simulations  (see 
Appendix  A).  Small  amplitude  ( v '  <  0.01  Uj)  random  velocity  fluctuations  were  added  to  this  in 
order  to  seed  the  turbulence.  The  consequence  of  not  adding  random  perturbations  was  a  prolonged 
region  of  laminar  flow  near  the  nozzle,  but  the  flow  downstream  was  not  particularly  sensitive  to  the 
nature  of  these  disturbances  provided  that  they  contained  a  range  of  frequencies  and  wavenumbers. 

Slot  jets  pictured  in  figure  0.1  were  used  to  excite  the  flow  in  a  manner  similar  to  the 
experiments  of  Parekh  et  al.7  Each  slot  extended  90°  around  the  jet  and  blew  normal  to  the 
shear  layers  just  downstream  of  the  nozzle.  The  techniques  for  including  these  actuators  in  the 
simulations  and  their  exact  specifications  are  outlined  with  the  numerical  method  in  Appendix  A. 
The  individual  slot  jets  blew  180°  out  of  phase  from  one  another  to  excite  a  flapping  mode  in 
the  jet  and  their  velocity  varied  sinusoidally  between  0  and  0.6Uj.  The  net  mass-flow  fraction  of 
the  actuators  was  Mact/Mjet  «  0.035.  Two  forced  jets  where  computed.  The  first  was  forced  at 
Strouhal  number  St  —  0.2,  which  had  been  found  experimentally  to  be  very  effective  at  spreading 
the  jet.7  The  other  case  was  forced  at  St  =  0.4,  which  was  the  most  amplified  frequency  of  the 
unforced  jet,  as  found  both  by  the  simulations  and  the  cited  experiment.14 


Results 


Visualizations 

Figure  0.2  shows  a  vorticity  magnitude  visualization  of  the  unforced  case.  The  initial  jet  shear  layers 
are  seen  to  be  laminar.  By  x  =  5r0  instability  waves  appear  which  develop  into  Kelvin-Helmholtz 
rollers  by  x  =  10ro.  Their  passing  frequency  is  St  =  0.4,  in  accord  with  the  experiments  of 
Stromberg  et  al.14  who  found  that  St  -  0.44  was  the  peak  Strouhal  number  in  the  early  development 
of  their  jet,  subject  to  ±10%  day-to-day  variation  in  their  facility.  At  the  instant  shown  the 
potential  core  extends  to  approximately  x  =  17 r0.  Near  the  end  of  the  potential  core  a  transition 
to  turbulence  occurs  which  is  corroborated  by  Reynolds  stress  statistics  that  will  be  reported 
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Figure  0.2:  Vorticity  magnitude  contours  for  the  unforced  case  in  an  x-r  plane.  There  are  40  evenly  spaced  contours 
with  maximum  ur0/Uj  —  8.43. 


x/r0 


Figure  0.3:  Actuator  flow  impinging  on  the  jet  shear  layer.  Every  other  mesh  point  in  x  and  every  sixth  mesh  point 
in  r  is  shown.  The  solid  rectangle  indicated  the  extent  of  the  actuator.  Light  contours  show  (p ~~ Poo) / PjUj  at  evenly 
spaced  intervals  (min  =  0.0333,  max  =  0.133). 


elsewhere.15 

Figure  0.3  shows  a  close-up  of  an  actuator  near  its  peak  blowing  condition.  (The  St  =  0.2 
case  is  shown,  but  the  St  =  0.4  case  is  similar.)  There  were  only  8  mesh  points  across  the  modeled 
actuator  in  the  streamwise  direction  and  those  near  the  center  of  the  actuator,  where  the  Gaussian 
axial  velocity  distribution  peaks  (see  Appendix  A ),  carry  most  of  the  momentum  flux.  A  region  of 
increased  pressure  just  below  the  actuator  is  also  shown  with  contours  in  the  figure.  The  initial 
effect  of  the  forcing  on  the  jet  is  seen  just  downstream  of  the  actuator  where  the  primary  jet  flow 
is  slightly  deflected.  The  bulk  of  the  fluid  exiting  the  actuator  appears  to  be  turned  downstream 
as  it  encounters  the  jet,  however  a  portion  of  it  is  also  turned  upstream,  giving  the  appearance  of 
a  stagnation  point  flow.  The  pressure  rise  in  this  region  is  also  reminiscent  of  a  stagnation  point 
flow. 

The  result  of  the  forcing  downstream  is  visualized  in  figures  0.4  and  0.5  for  the  St  =  0.2  and 
0.4  cases  respectively.  It  is  clear  for  the  St  —  0.2  case  that  the  actuators  excite  the  jet  significantly 
and  that  it  spreads  rapidly  in  the  9  =  7r/2  plane.  In  the  9  =  0  plane,  the  spreading  is  suppressed. 
A  single  mode  appears  to  dominate  the  early  development  of  the  jet  as  evidenced  by  the  large 
structures  visible  before  x  =  10 r0.  The  instantaneous  visualization  in  figure  0.4  shows  the  potential 
core  closing  at  approximately  x  =  9 r0.  Downstream  of  this  region,  the  vorticity  magnitude  field 
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Figure  0.4:  Vorticity  magnitude  contours  for  the  cased  forced  at  St  =  0.2.  Twenty  evenly  spaced  contours  between 
uiro/Uj  —  0.6  and  12.0  axe  shown. 


shows  an  eruption  of  small-scale  turbulence.  A  surprising  feature  of  the  visualization  is  the  apparent 
symmetry  near  x  -  7 r0  in  the  9  =  0  plane.  This  is  the  characteristic  feature  of  the  large-scale 
coherent  structures  seen  in  the  9  =  ir/l  visualization  as  they  intersect  the  0  =  0  plane. 

When  the  jet  is  forced  at  St  =  0.4  (figure  0.5),  the  large-scale  structures  are  smaller  but 
appear  earlier,  which  is  not  surprising  because  this  is  the  most  amplified  frequency  in  the  unforced 
jet.  However,  the  downstream  effect  of  the  forcing  is  now  quite  different.  The  structures  disappear 
or  are  obscured  by  small  scales  almost  immediately  and  the  jet  spreads  less  in  the  plane  of  the 
actuators  (9  =  n/2).  Though  the  shear  layers  appear  thicker  early  in  the  development,  they  slow 
their  spreading  and  merge  only  at  around  x  =  13r0,  beyond  where  they  merge  in  the  St  =  0.2  case. 
The  visualization  in  the  9  =  0  plane  is  similar  to  the  St  =  0.2  case,  with  similar  symmetries  at 
small  x  due  to  the  coherence  of  the  excited  structures,  but  a  longer  potential  core. 

Figure  0.6  shows  a  series  of  instantaneous  visualizations  of  the  jet  fluid  mixture  fraction  as 
it  adjusts  to  forcing  at  St  =  0.2.  The  time  series  starts  from  the  unforced  jet  and  the  actuators 
turn  on  at  t  =  0.  By  t  =  10ro/l7j,  there  is  clear  evidence  of  large  coherent  structures  distorting 
the  scalar  field.  As  expected,  these  travel  at  approximately  0.65 Uj,  appearing  at  x  «  10ro,  the 
horizontal  midpoint  of  the  region  shown,  at  t «  15r0/Uj.  The  scalar  field  takes  considerably  longer 
to  develop  further  downstream  where  decreasing  velocities  slow  advection.  For  computing  statistics, 
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Figure  0.5:  Vorticity  magnitude  contours  for  the  case  forced  at  St  =  0.4.  Twenty  evenly  spaced  contours  between 
< jJTo/Uj  =  0.6  and  12.0  are  shown. 
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t  =  0  r0/Uj 


t  =  10  r0/Uj  t  =  15  ro/Uj 


Figure  0.6:  Development  of  the  scalar  field  in  the  0  ==  7r/2  plane  once  St  —  0.2  forcing  is  turned  on  at  t  =  0.  Black  is 
pure  jet  fluid,  white  is  pure  ambient  fluid. 


the  forced  jets  were  assumed  to  be  fully  developed  only  after  t  =  80 r0/Uj.  The  simulations  were 
run  to  approximately  t  =  150ro /£/,-. 

Figure  0.7  shows  the  corresponding  set  of  images  for  the  jet  forced  at  St  =  0.4.  Again  we 
see  that  at  this  forcing  frequency,  large  structures  appear  closer  to  the  nozzle  than  in  the  St  =  0.2 
case.  The  mixed  regions  in  the  jet  shear  layers  thicken  rapidly,  but  a  tongue  of  pure  fluid  persists 
along  the  domain  centerline.  The  flapping  of  the  jet  column  is  also  less  pronounced. 

Mean  flow 

Mean  axial  velocity  (vx)  provides  a  more  quantitative  measure  of  the  effect  of  the  forcing  on  the 
jet  development.  Approximately  700  fields  spaced  in  time  by  A Tave  =  0.2 r0/Uj  were  averaged  to 
compute  the  mean.  Symmetries  were  exploited  to  increase  the  statistical  sample.  We  see  that 
the  jet  spreads  as  expected  (see  figure  0.8),  first  slowly  where  the  shear  layers  are  laminar,  and 
then  more  quickly  near  the  end  of  the  potential  core  where  the  flow  goes  through  transition.  If 
vx/Uj  =  0.9  is  used  to  designate  the  end  of  the  potential  core,  in  this  case  the  potential  core  closes 
at  x  —  17r0,  which  is  further  downstream  than  would  be  expected  for  high  Reynolds  number  jets, 
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t  =  10  ro/Uj  t  =  15  r0/Uj 


Figure  0.7:  Development  of  the  scalar  field  in  the  0  =  7r/2  plane  once  St  =  0.4  forcing  is  turned  on  at  t  =  0.  Black  is 
pure  jet  fluid,  white  is  pure  ambient  fluid. 
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Figure  0.8:  Mean  axial  velocity  ( vx/Uj )  for  the  unforced  jet. 


Figure  0.9:  Mean  axial  velocity  (vx/Uj)  for  St  =  0.2  forcing. 


because  the  shear  layers  are  initially  laminar  and  therefore  spread  slowly. 

Forcing  at  St  =  0.2  dramatically  alters  the  mean  flow  (figure  0.9).  In  the  plane  perpendicular 
to  the  action  of  the  actuators  (9  =  0)  the  jet  at  first  spreads  more  rapidly  than  the  unforced  case, 
but  then  this  is  reversed  starting  at  the  end  of  the  potential  core  ( x  =  8r0).  Only  near  x  =  20ro 
does  the  10  percent  velocity  contour  extend  to  the  same  radial  distance  as  at  x  =  8 r0.  The  contours 
have  an  unusual  thumb  shape  at  r  =  r0,  x  =  8r0  that  is  caused  by  the  large-scale  structures  seen  in 
figure  0.4.  Based  on  visualizations  ( e.g .  figure  0.4),  these  structures  first  intersect  the  9  —  0  plane 
near  r  =  0,  and  thus  they  bring  lower  velocity  fluid  into  that  region  before  the  region  near  r  =  r0. 
So  at  this  downstream  location  the  velocity  is  higher  near  r  =  r0  which  causes  the  appearance  of 
the  ‘thumb’.  In  the  9  =  n/2  plane,  the  jet  is  seen  to  spread  rapidly  starting  near  x  =  8r0  and 
this  continues  until  the  end  of  the  computational  domain  at  x  =  20 r0.  Parekh  et  al.7  showed  very 
similar  results  for  forcing  at  this  same  frequency. 

When  the  jet  is  forced  at  St  =  0.4  (figure  0.10),  the  mean  flow  is  markedly  different.  In 
the  0  =  0  plane,  the  jet  only  spreads  until  around  x  =  5.5 r0  before  spreading  is  reversed.  It  does 
this  significantly  closer  to  the  nozzle  than  in  the  St  —  0.2  case.  There  is  again  an  unusual  shape  to 
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Figure  0.10:  Mean  axial  velocity  ( vx/Uj )  for  St  =  0.4  forcing. 


the  contours  at  this  point,  but  at  this  forcing  frequency  it  is  less  pronounced.  Though  spreading 
in  the  &  =  0  plane  is  reversed  at  smaller  x  than  it  was  for  St  =  0.2  forcing,  the  potential  core  is 
now  longer,  extending  to  x  —  13.5r0.  In  the  6  =  7r/2  plane,  the  jet  spreads  rapidly  starting  at 
around  5.5r0,  but  spreading  slows  downstream  and  the  jet  does  not  grow  as  much  radially  as  in  the 
St  =  0.2  case  (figure  0.9). 


Unsteady  response  of  the  jet 

To  optimize  the  forcing,  by  either  an  open-  or  closed-loop  approach,  it  is  necessary  to  develop 
a  measure  of  its  effectiveness.  Naturally,  the  choice  of  a  definition  for  mixing  effectiveness  will 
depend  upon  the  specific  objective  of  the  application.  Here  we  will  consider  several  metrics  that 
are  of  potential  use  in  temperature  abatement  or  combustion  applications. 

Point  measurements  of  £ 

We  first  consider  point  measures  of  jet  fluid  mixture  fraction  (£  =  1  for  pure  jet  fluid  and  £  =  0 
for  pure  ambient  fluid)  on  the  jet  axis.  If  the  jet  were  hot,  the  concentration  of  jet  fluid  (mixture 
fraction)  would  closely  correspond  to  temperature.  Extensive  centerline  measurements  have  been 
made  in  forced  jets  and  have  been  used  to  estimate  the  mixing  enhancement  of  various  forcings.7 
Figure  0.11  shows  time  histories  of  (£)  on  the  jet  axis  (r  =  0)  and  at  x  =  5, 10, 15  and  20ro.  The 
angle  braces  ()  indicate  an  average  over  a  single  period  of  the  forcing:  Ta  =  9.Qr0/Uj  for  St  —  0.2 
and  Ta  —  A.9r0/Uj  for  St  =  0.4.  Despite  this  average,  there  is  still  considerable  oscillation  in  the 
measure  due  to  the  chaotic  nature  of  turbulence.  Averaging  for  longer  periods  would,  of  course, 
smooth  the  profiles,  but  for  closed-loop  control  applications  it  is  important  to  be  able  to  quickly 
measure  the  response  of  the  jet  to  changes  in  the  forcing.  Longer  averages  would  slow  the  response 
of  the  metric  to  the  changes  in  the  forcing. 

We  see  in  figure  0.11  that  before  forcing  is  initiated  at  t  =  0,  the  flow  is  pure  jet  fluid  at 
x  =  5 rQ  and  x  =  10ro.  At  x  =  15r0  there  is  slight  mixing  with  ambient  fluid  and  at  x  =  20, 
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Figure  0.11:  Time  histories  of  scalar  on  the  jet  axis: - x  =  5r0;  - x  =  10ro;  .  x  =  15r0;  and 

- x  =  20ro.  (a)  St  =  0.2;  (b)  St  =  0.4.  Each  curve  shows  a  time  average  over  one  forcing  period. 


the  mixture  fraction  hovers  around  its  long-time  (if  forcing  were  not  initiated)  mean  £  =  0.65. 
As  expected  from  observations  of  the  potential  core  length,  the  mixture  fraction  at  x  —  5 r0  is 
unaffected  by  the  forcing  at  either  Strouhal  number.  At  St  -  0.2  (figure  0.11  a),  the  mixture 
fraction  at  x  =  10ro  is  the  first  to  respond  to  the  forcing.  It  initially  decreases  to  (£)  =  0.5,  but 
rises  again  and  remains  near  the  mean  £  =  0.75  for  t  >  60ro/Uj.  The  period  averaged  values  at 
x  =  15r0  and  20ro  also  seem  to  overshoot  initially  before  they  settle  to  hover  around  their  apparent 
long-time  mean  values  of  £  =  0.6  and  £  =  0.3  respectively.  The  small  statistical  sample  size  makes 
these  values  and  the  point  where  they  are  reached  somewhat  imprecise.  Forcing  at  St  =  0.4  also 
causes  a  greater  reduction  in  centerline  mixture  fraction  initially  (figure  0.11  b).  The  curve  at 
x  =  10ro  initially  dips,  but  recovers  to  nearly  its  unforced,  pure  jet  fluid  level  by  t  =  70ro/Uj.  The 
curves  at  x  =  15r0  and  20ro  hover  around  their  mean  values  of  £  =  0.7  and  £  =  0.4.  It  is  unclear 
why  this  case  takes  longer  to  equilibrate  than  the  lower  frequency  forcing. 

The  initial  overreaction  of  the  jet  to  the  forcing  can  be  explained  qualitatively  with  a  simple 
model  where  the  large-scale  structures  are  assumed  to  be  linear  instability  waves.  Given  this  model, 
turbulent  structures  will  grow,  stabilize,  and  decay  as  the  layer  spreads.  If  the  jet  spreads  slowly, 
as  in  the.  unforced  case,  there  is  a  long  region  of  growth  before  decay.  So,  given  a  significant  initial 
forcing  amplitude,  the  structures  can  become  quite  intense  by  a  linear  mechanism.  However,  high 
amplitude  disturbances  will  also  increase  the  spreading  rate  of  the  layer  and  thereby  reduce  the 
distance  over  which  subsequent  disturbances  can  amplify.  Hence,  by  this  model  the  first  forced 
structure  “sees”  a  slowly  spreading  layer  and  thus  can  grow  more  than  subsequent  disturbances 
that  “see”  a  more  rapidly  spreading  layer,  which  explains  the  observed  overshoot. 

This  model  can  also  explain  the  greater  response  of  the  jet  to  St  =  0.2  forcing  than  forcing  at 
its  natural  frequency  St  =  0.4.  Instability  analysis  shows  that  for  thicker  shear  layers,  as  would  be 
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Figure  0.12:  Integrals  of  £"  with  n  =  1 - ;  n  =  2 - ;  n  —  4 .  ;  n  —  6  -  .  (a )  St  —  0.2  forcing;  (b) 

St  =  0.4  forcing. 

present  in  forced  jets,  the  most  amplified  instability  waves  will  have  a  longer  wavelength  and  lower 
frequency.16  Therefore,  it  is  not  surprising  that  St  =  0.2  is  more  successful.  Unfortunately,  linear 
stability  predictions  can  only  loosely  model  the  quantitative  behavior  of  turbulence.  An  accurate 
quantitative  prediction  using  linear  stability  analysis  is  therefore  not  likely  to  be  successful  in 
conjunction  with  the  high  amplitude  forcing  used  in  the  present  study. 

Volume  integrals  of  £n 

Volume  integrals  of  £n  can  also  provide  a  measure  of  mixing  effectiveness.  For  n  a  4  to  6  this  will 
provide  a  crude  model,  assuming  £  a  T,  of  the  infrared  signature  of  the  jet.  Time  histories  of 


where  Qc  is  the  physical  domain  (x  <  20 ra,r  <  10ro),  are  shown  in  figure  0.12.  No  time  averaging 
was  necessary  to  smooth  curves  in  this  case.  We  see  that  the  volume  integrals  of  £x  and  £2  both 
increase  when  the  jet  is  forced.  For  both  cases,  the  £J  curve  shows  that  there  is  60  percent  more  jet 
fluid  in  the  domain.  Despite  this  increase,  £4  and  £6  decrease  thus  indicating  improved  mixing.  The 
£6  curve  decreases  from  its  unforced  value  J0  to  J  =  0.7  J0  for  the  St  =  0.2  case  and  J  —  0.6  J0 
for  the  St  =  0.4  case.  It  is  somewhat  surprising  that  the  St  -  0.4  case  shows  better  mixing  by 
this  measure  because  the  opposite  was  predicted  based  upon  centerline  measurements  (figure  0.11). 
It  may  also  seem  counter  to  the  visualizations  and  mean  flow  data  which  show  greater  spreading 
for  St  =  0.2  forcing.  Because  (0.1)  depends  on  the  axial  dimension  of  the  computational  box,  it 
is  therefore  important  to  also  estimate  downstream  mixing  based  upon  the  available  data.  (Note 
that  the  finite  radial  dimension  does  not  affect  the  result  because  £  -»  0  by  r  =  10ro,  the  radial 
box  size.) 

We  can  make  such  an  extrapolation  by  computing  mixedness  as  a  function  of  downstream 
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Figure  0.13:  Planar  integrals  of  £6  from  (0.2):  - unforced; - St  =  0.2; .  St  —  0.4. 

distance.  Figure  0.13  shows  the  planar  contributions  to  J  as  a  function  of  x  for  t  >  80 r0/Uj: 
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It  now  becomes  clear  that  the  apparent  advantage  of  St  =  0.4  forcing  is  primarily  a  result  of  the 
jet’s  more  rapid  response  to  the  forcing.  At  the  outflow  boundary  we  see  that  the  St  =  0.2  forced 
jet  is  actually  mixed  better  (by  this  measure)  and  Mx  has  a  steeper  slope  (rate  of  mixing)  than 
the  jet  forced  at  St  =  0.4.  Though  it  is  not  possible  to  make  firm  judgments  about  the  subsequent 
downstream  mixing,  based  on  the  level  and  slope  of  Mx  at  x  =  20 r0  it  appears  that  the  St  =  0.2 
case  might  be  better  if  more  downstream  fluid  could  be  included  in  (0.1).  Both  forced  cases  are 
clearly  better  than  the  unforced  case  also  shown  in  the  figure. 

The  asymmetry  of  the  jets  is  seen  in  figure  0.14  which  shows  planar  integrals  of  £6  at  constant 

0 , 

/20ro  rl0ro 

/  £ fdrdx .  (0.3) 

-20 r0  Jo 

For  both  forced  cases,  this  metric  peaks  in  the  9  =  7r/2  plane,  the  plane  of  the  actuators.  For 
both  Strouhal  numbers,  this  peak  is  roughly  1.75  times  as  high  as  the  minimum  values.  Somewhat 
surprisingly  these  minima  do  not  occur  at  9  =  0  (see  figure  0.14)  which  is  perpendicular  to  the 
actuators.  Mg  is  20  percent  higher  at  9  =  0  than  at  its  minima.  Nearly  everywhere  Mg  is 
suppressed  below  its  value  in  the  unforced  case,  more  so  in  the  St  =  0.4  case.  It  should  be  noted 
that  the  details  of  these  results  also  depend  upon  the  box  size.  Because  the  forced  jets  are  still 
highly  asymmetric  at  the  outflow  boundary,  the  numerical  differences  between  the  peaks  in  valley 
would  increase  if  more  downstream  data  were  available. 

Streamwise  mass  flux 

The  net  entrainment  of  the  jet  can  be  studied  by  computing  the  streamwise  mass  flux.  Because  vx 
is  negligible  at  r  =  10ro,  this  is  equivalent  to 

/*2  7T  plO  r0 

F(x)=  /  pv^rdrd9,  (0.4) 

Jo  Jo 

which  is  plotted  in  figure  0.15.  The  mass  flux  in  the  unforced  case  grows  slowly  at  first,  where 
the  shear  layers  are  laminar  and  disturbances  are  small,  and  then  grows  more  rapidly  starting  at 
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Figure  0.14:  Planar  integrals  of  f6  from  (0.3):  - unforced; - St  =  0.2; .  St  =  0.4. 


around  x  =  14 r0  where  the  potential  core  closes.  The  mass  fluxes  in  the  forced  cases  both  grow 
rapidly  from  the  start  and  are  over  twice  at  high  by  x  —  20 r0.  Forcing  at  St  =  0.2  is  only  mildly 
more  effective  at  increasing  the  mass  flux  than  forcing  at  St  —  0.4.  It  is  seen  that  the  fluxes  are 
slightly  different  for  the  three  cases  at  x  =  0.  This  is  because  of  the  fluid  added  by  the  actuators 
and  entrainment  caused  directly  by  their  action. 

Scalar  dissipation 

The  rate  of  scalar  dissipation  can  also  provide  an  important  measure  of  mixing  which  is  particularly 
relevant  in  combustion  applications.  In  figure  0.16,  we  consider  planar  integrals  of  |V£|  as  a  function 
of  downstream  distance, 

r2ir  rlQr0 

V{x)=  /  |V£|  rdrdd.  (0.5) 

Jo  Jo 

V{x)  for  the  St  -  0.4  case  starts  to  rise  closer  to  the  nozzle,  but  the  curve  for  St  =  0.2  forcing 
follows  only  2r0  further  downstream  and  becomes  50  percent  greater  by  the  right  side  of  the  domain. 
The  total  dissipation  in  the  computational  domain  is  clearly  larger  for  St  =  0.2  forcing.  Based  on 
the  data  at  x  —  20 r0,  this  trend  is  likely  to  continue  downstream. 


Conclusion 

Simulations  of  jets  forced  with  high  amplitude  actuation  reproduced  experimental  observations  of 
similarly  forced  jets.  Visualizations  showed  dramatic  effect  of  this  forcing  on  the  jet.  Forcing  at 
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Figure  0.16:  Planar  integrals  of  |V£|  from  (0.5):  - unforced; - St  =  0.2; .  St  —  0.4. 

St  =  0.2  excited  the  jet  column  into  a  distinct  flapping  mode.  When  forced  at  St  =  0.4,  the  most 
amplified  frequency  in  the  unforced  jet,  the  large  structures  appeared  closer  to  the  nozzle.  However, 
despite  the  rapid  initial  spreading  of  this  jet,  forcing  at  St  =  0.4  was  not  as  effective  at  spreading 
the  jet  and  mixing  it  (by  most  measures)  with  the  ambient  flow  downstream.  Mean  axial  velocities 
showed  that  the  jet  becomes  highly  non-axisymmetric  in  both  forced  cases  and  that  the  potential 
core  length  was  reduced  more  by  St  =  0.2  than  St  =  0.4  forcing. 

Mixing  was  quantified  by  several  different  metrics.  Point  measurements  of  scalar  concentra¬ 
tion  on  the  jet  axis  showed  that  forcing  at  St  =  0.2  was  more  effective  at  reducing  centerline  scalar 
concentration.  However,  volume  integrals  of  £”  over  the  computational  domain  (for  n  =  4  or  6) 
were  smaller  when  the  jet  was  forced  with  St  =  0.4.  This  was  primarily  due  to  the  faster  response 
of  the  jet  to  forcing  at  this  frequency.  Consideration  of  planar  integrals  of  £6  and  |V£|  as  a  function 
of  x  suggested  that  forcing  at  St  =  0.2  would  be  more  effective  for  mixing  further  downstream. 
Forcing  at  both  Strouhal  numbers  increased  streamwise  mass  flux  considerably  over  the  unforced 
case,  with  St  —  0.2  forcing  performing  marginally  better  in  this  regard  than  St  =  0.4  forcing. 
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Appendix  C 

Optimization  of  Turbulent  Jet  Mixing 


Summary 

We  introduce  an  approach  for  controlling  jet  mixing  that  combines  direct  numerical  simulation  of 
an  incompressible  jet  flow  with  stochastic  optimization  procedures.  The  jet  is  excited  with  helical 
and  combined  helical  and  axial  actuations  at  the  orifice.  An  objective  function  that  measures 
the  spreading  of  the  jet  evaluates  the  performance  of  the  actuation  parameters.  The  optimization 
procedure  searches  for  the  best  actuation  by  automatically  varying  the  parameters  and  calculating 
their  objective  function  value.  Solutions  that  lead  to  a  pronounced  spreading  of  the  jet  are  found 
within  reasonable  time,  although  the  evaluation  of  the  objective  function,  the  DNS  or  LES  of  the  jet, 
is  expensive.  For  a  jet  flow  at  low  Reynolds  number  the  performance  of  different  search  algorithms 
(simulated  annealing  and  evolution  strategies)  is  evaluated.  We  compare  various  objective  functions 
based  on  radial  velocity  and  the  concentration  of  a  passive  scalar,  including  functions  that  penalize 
actuation  with  high  amplitudes.  We  find  that  a  combined  axial  and  helical  actuation  is  much  more 
efficient  with  respect  to  jet  mixing  than  a  helical  actuation  alone. 
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Introduction 

The  control  of  turbulent  jet  flows  has  applications  in  various  fields  such  as  combustion,  aerodynamic 
noise  and  jet  propulsion.  In  combustion  processes  it  is  important  to  enhance  the  turbulent  mixing  of 
the  chemical  species  to  make  the  process  more  efficient  and  reduce  the  concentration  of  pollutants. 
Acoustic  emissions  by  aircraft  jet  engines  can  be  reduced  by  controlling  flow  unsteadiness  that 
produces  noise.  One  aim  of  enhanced  mixing  in  jet  propulsion  applications  is  to  decrease  the  plume 
temperature  and  suppress  infrared  radiation. 

The  mixing  rate  of  a  turbulent  jet  can  be  significantly  altered  by  applying  a  suitable  excitation 
at  the  jet  orifice.  Since  the  external  forcing  interacts  with  the  natural  modes  of  the  jet  in  a  nonlinear 
way,  it  is  not  possible  to  predict  which  kind  of  actuation  is  optimal  to  increase  mixing.  Various 
experiments  have  been  carried  out  that  study  the  reaction  of  jets  to  the  nozzle  geometry  (passive 
control  [1])  and  to  external  forcing  (active  control  [2,  3,  4]).  Different  types  of  actuators  such  as 
piezoelectric  devices  and  synthetic  jets,  have  been  tested  [5,  6].  It  has  been  shown  that  a  large 
spreading  of  the  jet  can  be  achieved  with  a  small  mass  flow  actuation  if  suitable  frequencies  are 
chosen  [6]. 

Numerical  simulations  of  compressible  and  incompressible  jet  flows  have  been  carried  out  that 
confirm  many  observations  made  in  experiments  with  periodically  forced  jets  [7,  8,  9,  10].  So  far 
it  is  difficult  to  carry  out  a  systematic  search  for  the  optimal  forcing  because  the  simulations  are 
very  computationally  intensive.  Evolution  strategies  have  been  first  applied  to  the  optimization  of 
actuation  parameters  in  jet  mixing  by  Koumoutsakos  et  al.  [11].  This  work  showed  that  evolu¬ 
tion  strategies  are  capable  of  finding  suitable  actuations  for  a  vortex  model  and  direct  numerical 
simulations  of  compressible  jets. 

In  this  paper  we  perform  a  systematic  search  for  helical  and  combined  axial  and  helical  forcing  of 
a  jet  that  maximize  mixing,  as  determined  by  various  cost  functions.  We  combine  automatic  search 
strategies  with  direct  numerical  simulation  (DNS)  and  large  eddy  simulation  (LES)  of  a  round  jet. 
Varying  the  actuation  at  the  orifice,  we  search  for  the  forcing  that  maximizes  the  spreading  of  the 
jet.  Large  spreading  corresponds  to  efficient  mixing  of  a  passive  scalar,  transported  by  the  flow,  with 
the  surrounding  air.  The  phase  space  of  possible  actuation  parameters  is  searched  automatically 
by  a  stochastic  optimization  strategy.  While  classical  approaches  like  gradient  methods  converge 
rapidly,  there  is  a  risk  of  premature  convergence  to  a  local  optimum.  Stochastic  methods  avoid  this 
by  allowing  steps  in  seemingly  unfavorable  directions  with  a  certain  probability.  They  are  therefore 
useful  for  the  optimization  of  multimodal  functions.  Among  the  stochastic  methods  proposed  in 
the  literature  [12],  simulated  annealing  and  evolution  strategies  have  been  used  for  a  large  number 
of  engineering  optimization  problems  [13].  Simulated  annealing  [14,  15],  also  known  as  Monte  Carlo 
annealing,  is  a  serial  procedure  which  is  very  efficient  and  straight  forward.  Evolution  programs 
is  the  name  for  a  group  of  strategies  which  make  use  of  the  principles  of  evolution  -  reproduction, 
mutation  and  selection  -  to  find  an  optimal  solution  [16,  17,  18,  19],  They  are  easy  to  implement, 
efficient  and  inherently  parallel.  Examples  for  the  application  of  evolution  strategies  to  problems 
in  fluid  mechanics  and  first  results  for  the  optimization  of  jet  mixing  have  been  presented  in  [20]. 
In  this  paper,  we  will  give  a  detailed  description  of  the  optimization  procedure  and  compare  the 
performance  of  different  search  strategies  for  the  optimization  of  jet  actuation.  A  fast  convergence 
of  the  procedure  is  indispensable  since  the  numerical  simulation,  which  evaluates  a  given  actuation, 
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is  very  computationally  expensive.  In  order  to  compare  the  different  strategies  we  have  chosen  a 
simple  test  case,  a  jet  at  low  Reynolds  number,  for  which  the  optimal  solution  is  known. 

Direct  Numerical  Simulation  of  the  Jet 

The  numerical  solver  simulates  a  free  round  jet  issuing  from  a  circular  orifice  of  diameter  D  in  a 
solid  wall.  The  incompressible  Navier-Stokes  equations  and  the  transport  equation  for  a  passive 
scalar  are  solved  in  a  spherical  coordinate  system,  with  (r,  0,  <f>)  denoting  the  radial,  tangential  and 
azimuthal  directions  (Fig.  1).  Spatial  derivatives  are  calculated  on  a  staggered  spherical  grid.  The 
time  integration  is  carried  out  with  a  second  order  Adams-Bashforth  method.  More  details  of  the 
numerical  scheme  can  be  found  in  Boersma  [21].  The  computational  domain  shown  in  Fig.  1  results 
from  the  intersection  between  the  shell  defined  by  the  surfaces  r  =  5D  and  r  =  15D  and  the  cone 
starting  from  the  center  of  the  sphere  with  an  opening  angle  of  36°.  This  geometry  covers  a  domain 
with  a  streamwise  extent  of  10D  and  a  spanwise  diameter  of  4 D  for  the  inflow  section  and  10 D  for 
the  outflow  section.  Such  a  discretization  is  able  to  follow  the  streamwise  spreading  of  the  jet  and 
to  allow  a  well-balanced  resolution  of  the  flow  field  with  a  reasonable  number  of  grid  points. 

At  the  inflow  section,  the  following  mean  streamwise  velocity  profile  is  imposed 

Vzo(rc)  =  +tanh[0.25D/Q0(D/(4rc)  -4rc/£>)]),  (1) 

where  ho  is  the  centerline  velocity  and  rc  is  the  radius  in  a  cylindrical  system  ( z,rc,</> ).  The  initial 
momentum  thickness  was  0o  =  -D/60  in  our  simulations. 

At  the  lateral  boundary,  the  total  normal  stress  is  set  to  zero,  which  allows  fluid  exchange  across 
the  boundary.  This  condition  properly  simulates  the  entrainment  of  ambient  fluid  in  the  spreading 
jet  flow.  A  so-called  convective  boundary  condition  (Orlanski  [22])  is  used  to  evacuate  the  vortex 
structures  through  the  downstream  boundary.  This  condition  is  numerically  stable  but  physically 
not  very  realistic  in  elliptic  flows  like  the  jets  described  here.  However,  the  convective  nature  of 
the  homogeneous  jet  flow  ([23])  allows  the  eventual  removal  of  spurious  disturbances. 

For  the  DNS  of  a  jet  with  Reynolds  number  Re  —  1500,  based  on  orifice  velocity  Vo  and  diameter 
D,  the  spherical  grid  consists  of  192  x  128  x  96  points  in  the  radial,  tangential  and  azimuthal 
directions  respectively.  Typical  CPU  times  for  the  calculation  of  a  fully  developed  jet  are  1200  node 
hours  on  an  Origin  2000,  using  the  Message  Passing  Interface  (MPI).  The  optimization  requires 
the  simulation  of  approximately  150  jets  for  different  actuation  parameters.  We  will  show  in  in  this 
paper  how  we  have  reduced  the  CPU  time  during  the  optimization  process. 

In  unforced  jets,  large  coherent  structures  are  observed  that  are  related  to  the  instability  modes 
of  the  jet.  The  dominant  modes  are  the  axisymmetric  or  varicose  mode  and  helical  modes.  The 
axisymmetric  mode  causes  the  shear  layer  to  roll  up  into  vortex  rings.  By  applying  axial  forcing 
to  the  shear  layer,  the  frequency  of  vortex  ring  generation  and  the  pairing  of  the  vortex  rings  may 
be  altered.  The  initial  shear  layer  is  able  to  amplify  a  large  range  of  frequencies.  The  frequency  / 
which  leads  to  the  maximum  amplification  of  the  initial  shear  layer  is  called  the  natural  frequency. 
It  can  be  obtained  by  linear  spatial  instability  analysis  (Michalke  [24],  Ho  and  Huerre  [25]).  For 
an  axisymmetric  jet  with  the  initial  velocity  profile  given  in  Eqn.  (1)  it  is  StQ0  —  fQo/Vo  —  0.018. 
The  frequency  fp  of  the  axial  perturbation  that  produces  the  largest  total  amplification  is  called 
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the  preferred  mode  of  the  jet.  It  corresponds  to  the  frequency  of  vortices  at  the  end  of  the  potential 
core.  This  frequency  has  been  determined  by  Crow  and  Champagne  to  be  Stp  «  0.3  [261.  However, 
other  studies  have  found  the  preferred  Strouhal  number  to  vary  between  0.25  and  0.5  [27,  28!. 
Mankbadi  [28]  observed  for  a  round  jet  under  axisymmetric  forcing  that  mixing  is  enhanced  if  the 
forcing  Strouhal  number  is  equal  to  about  twice  the  jet’s  preferred  mode.  It  was  also  found  that 
at  high  Strouhal  numbers  the  momentum  thickness  is  reduced  along  the  jet. 

Mixing  can  be  increased  significantly  if  axial  and  helical  forcing  are  combined  [2,  3,  6].  In  the 
following  we  will  denote  the  Strouhal  numbers  of  axial  and  helical  forcing  with  Sta  and  Sth .  Large 
spreading  angles  (up  to  80°)  have  been  observed  for  certain  ratios  (3  of  the  Strouhal  numbers  Sta 
and  Sth  [4].  In  particular,  the  jet  splits  into  two  branches  for  0  =  2  (bifurcating  jet)  and  into 
three  branches  for  0  =  3,  For  non-integer  ratios  1.6  <  /?  <  3.2,  vortex  rings  are  shed  in  different 
directions  (blooming  jet) .  It  has  been  discussed  by  Parekh  et  al.  [4]  that  these  flow  patterns  appear 
for  both  low  and  high  Reynolds  numbers. 

In  our  simulation  the  actuation  of  the  shear  layer  was  achieved  by  superposing  periodic  distur¬ 
bances  on  the  initial  velocity  profile.  We  have  used  two  kinds  of  actuation.  For  the  first,  the  total 
inflow  velocity  in  the  z  direction  of  a  cylindrical  coordinate  system  (z,rc,<p)  is 


Vz{rc,(p,t)  =  Vz0  (r) 


/  y0  \  2  rc 

1  +  Ahsin  (2 TrSth—  t)  cos(ip)  — 


(2) 


Ah  is  the  amplitude  of  the  actuation  which  is  phase- locked  in  the  plane  <p  —  0.  It  corresponds  to  the 
superposition  of  two  counter-rotating  helical  modes  of  the  same  frequency.  It  has  been  observed 
in  simulations  and  experiments  that  this  type  of  actuation  causes  the  jet  to  perform  a  flapping 
motion  in  the  plane  of  actuation.  The  second  type  of  actuation  is  a  superposition  of  axial  arid 
helical  modes 


Vz(rc,tp,t)  =  Vz0  (r) 


1  +  Aasin  (2n StJjj  tj  +  Ahsin  (ivSth^t  +  a'j  cos(tp) 


2  rc 

D 


(3) 


Here  Aa  and  Ah  are  the  amplitudes  of  the  axial  and  helical  mode  and  Sta,Sth  are  the  respective 
Strouhal  numbers.  We  have  again  chosen  the  helical  part  of  the  actuation  to  be  phase  locked  in 
the  plane  <p  =  0.  The  angle  a  determines  the  plane  of  spreading.  The  actuation  Eqn.  (3)  with 
(3  =  Sta/ Sth  =  2  has  been  used  before  to  model  bifurcating  jets  [21].  To  simplify  the  notation,  we 
will  refer  to  jets  obtained  with  the  dual-frequency  actuation  Eqn.  (3)  as  bifurcating  jets  if  f}  takes 
on  a  value  that  is  approximately  two.  In  the  following  sections  we  will  investigate  which  parameter 
vector  x  =  (Sta,  Sth,  An  Ah)  maximizes  the  spreading  of  the  jet. 


Optimization  Strategies 

The  optimization  problem  is  described  by  a  vector  of  parameters  that  are  varied  during  the  proce¬ 
dure,  and  an  objective  function  that  evaluates  the  performance  of  the  parameters.  In  the  following, 
we  will  use  the  terminology  of  evolution  programs.  We  will  refer  to  the  objective  function  as  fitness 
function,  to  its  value  after  evaluation  as  fitness  value,  and  to  the  surface  that  describes  the  fitness 
value  as  a  function  of  the  parameters  as  fitness  landscape.  We  will  call  the  sequence  of  parameter 
vectors  that  are  evaluated  during  the  optimization  procedure  a  search  trajectory.  For  our  opti¬ 
mization  problem,  the  parameter  vector  consists  of  the  Strouhal  numbers  and  amplitudes  of  the 
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actuation.  The  objective  function  measures  the  spreading  of  the  jet  by  performing  a  direct  numer¬ 
ical  simulation  of  the  jet  flow.  The  fitness  landscape  hence  describes  how  jet  spreading  depends  on 
the  frequencies  and  amplitudes  of  the  actuation. 

For  problems  with  unimodal  fitness  landscapes,  classical  optimization  strategies  like  the  gradient 
method  are  very  efficient.  For  problems  with  multimodal  fitness  landscape  however  there  is  a 
risk  for  deterministic  methods  to  converge  to  a  local  optimum  (so-called  premature  convergence). 
Stochastic  methods  do  not  always  follow  the  gradient  during  the  optimization  process,  but  instead 
allow  large  search  steps  as  well  as  steps  in  directions  that  may  have  a  worse  fitness  value  with  a 
certain  probability.  This  enables  the  search  trajectory  to  leave  a  local  optimum.  On  the  other 
hand,  these  methods  usually  require  large  numbers  of  iterations  when  compared  with  deterministic 
gradient-based  methods.  The  dynamic  behavior  of  the  jet  is  determined  by  the  nonlinear  interaction 
of  different  modes.  This  is  likely  to  cause  a  complicated  dependence  of  the  jet  development  on  the 
actuation  parameters.  We  therefore  expect  the  objective  function  to  be  multimodal,  which  suggests 
that  stochastic  procedures  are  the  method  of  choice  for  optimization  of  jet  mixing. 

Our  approach  of  optimizing  jet  mixing  implies  that  for  each  search  step  a  direct  numerical 
simulation  of  the  jet  flow  has  to  be  performed,  which  is  very  expensive.  Even  when  various  ap¬ 
proximations  to  reduce  the  CPU  time  and  parallel  processing  of  the  jet  code  are  used,  the  total 
number  of  jet  flow  calculations  that  can  be  achieved  is  of  the  order  100-200.  For  most  optimization 
problems,  typical  numbers  of  objective  function  evaluations  are  of  the  order  of  103  —  10°.  A  rapid 
convergence  of  the  search  is  therefore  indispensable. 

We  will  compare  the  following  optimization  procedures:  a  simulated  annealing  algorithm,  im¬ 
plemented  by  Goffe  et  al.  [29],  a  single- membered  and  a  multi-membered  evolution  strategy  with 
Covariance  Matrix  Adaption  of  the  steplength  [30],  both  implemented  by  Koumoutsakos  and  Muller 
[20].  The  single- membered  evolution  strategy  corresponds  to  a  simple  stochastic  search,  which  only 
accepts  steps  that  lead  to  an  improved  fitness  value,  while  the  simulated  annealing  algorithm  allows 
uphill  movement  to  some  extent  (for  a  minimization  problem).  Both  methods  are  serial  and  use 
step  length  adjustment  based  on  the  fraction  of  successful  steps.  The  multi-membered  evolution 
strategy  operates  on  larger  populations  of  individuals  and  evaluates  multiple  trajectories  in  parallel. 
It  includes  a  refined  step  length  adaption  scheme  which  makes  use  of  the  whole  search  path  rather 
than  just  the  number  of  successful  trials. 

Evolution  strategies 

Evolutionary  computation  techniques  employ  operations  inspired  by  biological  principles  such  as 
a  reproduction  cycle,  natural  selection  and  diversity  by  variation,  to  identify  the  optimal  solution. 
Examples  are  genetic  algorithms  [16],  evolutionary  programming  [18]  and  evolution  strategies  (ES) 
[19].  It  has  been  shown  by  Koumoutsakos  et  al.  [11]  that  evolution  strategies  are  well  suited  for 
the  optimization  of  jet  actuation. 

Evolution  strategies  operate  on  populations  of  individuals  which  represent  possible  solutions 
of  the  given  problem.  Each  individual  consists  of  a  vector  x  containing  n  parameters  and  an 
associated  vector  s  containing  n  mutation  steplength.  The  number  of  parents  in  each  generation  is 
p,  the  number  of  offspring  is  A,  with  A  >  p.  Starting  from  an  initial  generation  of  p  parents,  the 
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offspring  are  generated  by  mutation 


(4) 


x*+1=4  +  spV(0,l). 

Here,  the  indices  p  and  o  denote  parent  and  offspring  and  i  is  the  number  of  the'  generation. 
jV(0, 1)  is  a  normal  distributions  with  zero  average  and  unit  variance.  We  impose  certain  limits  to 
the  n  parameters.  The  amplitudes  must  be  positive  and  the  Strouhal  numbers  should  be  within  a 
certain  interval  (see  below).  In  our  optimization,  a  mutation  step  that  leads  to  a  value  outside  the 
allowed  range  is  considered  unsuccessful  and  repeated  until  allowed  values  are  obtained.  If  several 
repetitions  of  the  mutation  step  fail,  the  parameters  are  defined  randomly  within  the  given  limits. 

After  A  offspring  have  been  generated,  their  fitness  values  are  determined  and  the  best  p  in¬ 
dividuals  are  selected.  The  steplengths  of  the  mutation  are  decreased  if  the  offspring  are  better 
than  the  parents,  otherwise  they  are  increased.  The  p  best  individuals,  i.e.  the  vectors  of  the  best 
parameters  and  the  corresponding  mutation  steplengths  are  used  as  parents  for  the  next  generation. 
The  best  individuals  are  chosen  among  the  offspring  {{p,  A)  strategy)  or  among  the  p  +  A  parents 
and  offspring  ((p  +  A)  strategy)  [31].  The  p  best  individuals,  i.e.  the  vectors  of  the  best  param¬ 
eters,  may  also  inherit  their  associated  mutation  steplength.  The  simplest  and  earliest  version  of 
the  evolution  strategy  [19]  describes  a  population  consisting  of  just  one  parent  and  one  offspring 
and  is  referred  to  as  (1  +  1)  strategy.  It  corresponds  to  a  stochastic  search  which  accepts  a  new 
solution  only  if  it  has  a  better  objective  function  value  than  the  previous.  It  uses  a  step  length 
adjustment  based  on  the  ratio  of  the  number  of  successful  to  the  total  numbers  of  trials.  The  use  of 
several  search  trajectories  allows  the  simultaneous  evaluation  of  different  optima  and  therefore  can 
speed  up  the  search  significantly.  Hansen  et  al.  [30,  32]  have  shown  that  higher  convergence  rates 
can  be  achieved  using  the  Covariance  Matrix  Adaptation  scheme  (CMA).  With  this  method,  the 
step  sizes  are  adapted  by  considering  the  length  of  successful  mutation  steps  not  only  of  the  last 
but  rather  of  a  number  of  previous  generations.  The  CMA  of  step  sizes  makes  use  of  correlations 
between  consecutive  steps.  If  successive  mutation  steps  are  positively  correlated,  they  should  be 
replaced  by  one  larger  step.  If  the  correlation  is  negative,  the  steps  partly  cancel  and  should  be 
replaced  by  a  shorter  step.  The  procedure  is  thus  most  efficient  if  there  is  no  correlation  between 
consecutive  steps,  which  means  that  they  are  orthogonal.  In  order  to  achieve  this,  an  evolution  path 
is  calculated  as  the  sum  of  weighted  steplength.  From  the  evolution  path  a  covariance  matrix  is 
obtained,  which  is  used  to  turn  the  mutation  vector  in  the  most  promising  direction.  The  adaption 
of  the  evolution  path  takes  place  on  a  rather  large  time  scale.  In  order  to  allow  fast  changes  of 
the  step  lengths,  a  global  step  length  is  introduced  which  changes  on  a  short  time  scale.  A  further 
speed-up  is  achieved  by  combining  the  CMA-ES  with  an  intermediate  recombination  that  averages 
the  variable  vector  elements  of  some  of  the  parents  [30].  Details  of  the  CMA-ES  strategy  can  be 
found  in  [30].  The  search  is  terminated  if  either  the  improvement  made  during  the  search  gets 
smaller  than  a  certain  value  or  if  a  maximum  number  of  function  evaluations  is  reached. 

Simulated  Annealing 

Simulated  annealing  is  based  on  the  analogy  between  the  annealing  of  solids,  i.  e.  the  way  in  which 
a  solid  cools  and  freezes  into  a  crystalline  structure  with  minimal  energy,  and  the  problem  of  solving 
large  optimization  problems.  The  algorithm,  first  proposed  by  Kirkpatrick  et  al.  [14]  is  easy  to 
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implement  and  very  efficient.  It  employs  a  random  search  which  does  not  only  accept  changes  that 
improve  the  objective  function  /,  such  that  A /  =  /(xn+1)  -  /(xn)  <  0  in  a  minimization  problem. 
It  also  accepts  an  increase  A/  >  0  with  probability 


p  —  exp 


(5) 


where  T  is  the  ’system  temperature’  which  is  high  in  the  beginning  of  the  annealing  process  and 
approaches  zero  towards  its  end.  Eqn.  (5),  known  as  the  Metropolis  criterion  [33],  determines  the 
amount  of  uphill  movement  that  is  allowed.  The  optimization  procedure  consists  of  a  sequence  of 
’temperatures’  T.  At  each  T,  a  random  search  with  a  fixed  number  of  search  steps  (a  Markov  chain 
of  certain  length)  is  generated,  where  new  parameter  vectors  are  obtained  similar  to  Eqn.  (4).  If  no 
better  solution  has  been  found  at  the  end  of  the  chain,  the  start  parameters  of  the  search  are  kept 
as  initial  points  for  the  search  at  the  next  T.  The  parameter  T  ensures  that  at  the  beginning  of  the 
search  steps  are  possible  that  lead  to  a  worse  objective  function  value,  which  enables  the  trajectory 
to  leave  a  local  optimum.  At  the  end  of  the  search,  when  the  vicinity  of  the  global  optimum 
is  reached,  the  probability  of  accepting  a  worse  solution  is  very  low  to  ensure  rapid  convergence 
to  the  optimum.  As  rule  for  decrementing  the  temperature,  the  so-called  annealing  schedule,  we 
have  chosen  an  exponential  cooling  scheme  [14].  The  initial  temperature  was  set  such  that  the 
average  probability  for  accepting  a  worse  parameter  vector  is  about  0.7  [34].  The  steplength  is 
again  adjusted  according  to  the  successful  number  of  trials. 


Jet  Optimization 

In  this  section,  we  describe  the  optimization  procedure  and  compare  the  performance  of  the  different 
algorithms.  We  have  to  choose  an  objective  function  which  provides  an  effective  measure  of  the 
spreading  of  the  jet  and  which  is  sufficiently  sensitive  to  changes  in  the  actuation  parameters.  It  is 
important  to  ensure  that  the  evaluation  of  the  objective  function  is  possible  within  reasonable  CPU 
times.  We  introduce  two  approximations  that  lead  to  a  reduction  of  the  CPU  time:  The  evaluation 
of  the  objective  function  at  an  early  stage  of  the  jet  simulation  and  the  use  of  a  coarse  grid  for 
these  calculations.  We  investigate  a  test  case,  a  jet  at  low  Reynolds  number  Re  =  1500,  for  which 
the  phase  space  can  be  examined.  First  results  for  jets  at  higher  Reynolds  number,  Re  =  6000,  are 
presented. 

Optimization  Parameters 

The  actuation  of  the  jet  is  described  by  the  parameter  vector  x.  If  the  single-frequency  actua¬ 
tion  Eqn.  (2)  is  chosen,  the  vector  is  two-dimensional  and  consists  of  the  Strouhal  number  and 
amplitude  x  =  (St,  Ah)-  For  the  dual-frequency  excitation  Eqn.  (3)  the  vector  consists  of  two 
Strouhal  numbers  and  two  amplitudes  x  =  ( Sta,Sth,Aa,Ah ),  or  x  =  (Sta,P,Aa,Ah)  if  the  ratio 
of  frequencies  fi  =  Sta/Sth  is  used.  The  limits  of  the  parameters  have  been  chosen  as  follows.  In 
laboratory  experiments,  bifurcating  jets  have  been  observed  for  0.4  <  Sta  <  0.65,  (3  =  2  [4].  In 
experiments  by  Parekh  et  al.  [6]  that  used  an  actuation  similar  to  Eqn.  (2)  it  was  found  that  the 
most  pronounced  flapping  of  the  jet  appears  for  approximately  Sth  =  0.2.  A  large  amplification  of 
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signals  was  also  observed  for  St  ss  0.27  and  St  «  0.4.  Numerical  simulations  by  Freund  and  Moin 
[8]  have  shown  that  flapping  jets  appear  for  Sth  =  0.2  and  Sth  =  0.4,  and  that  the  lower  Strouhal 
number  is  more  efficient  for  jet  mixing.  This  Strouhal  number  has  also  been  found  in  the  first 
attempts  to  find  optimal  actuation  of  a  compressible  jet  with  an  evolution  strategy  [11].  We  have 
varied  the  Strouhal  numbers  in  a  somewhat  wider  range  0.1  <  Sth  <  0-55  for  the  single- frequency, 
and  0.1  <  Sta  <  1.2  for  the  dual-frequency  excitation.  For  the  latter  case,  the  ratio  of  Strouhal 
numbers  was  1.6  <  (3  <  3.2. 

In  laboratory  experiments  jets  are  perturbed  by  acoustic  forcing,  in  recent  experiments  and 
applications  by  piezoelectric  devices  or  by  suction  and  blowing  of  fluid  near  the  jet  orifice  (syn¬ 
thetic  jets).  While  acoustic  actuation  with  high  amplitude  cannot  be  easily  achieved,  blowing  and 
suction  with  a  large  mass  flow  is  not  practical  and  should  therefore  be  avoided.  In  order  to  match 
this  situation  we  need  to  search  for  an  actuation  that  leads  to  a  large  spreading  of  the  jet  while 
maintaining  a  sufficiently  low  mass  flow.  In  part  of  our  optimization  runs  we  have  included  this 
constraint  by  imposing  an  upper  limit  Ah  <  0.06  and  by  introducing  a  suitable  penalty  term  in  the 
objective  function. 

Choice  of  the  Objective  Function 

The  objective  function  measuring  the  spreading  of  the  jet  is  based  on  quantities  available  by  the 
DNS  simulation:  the  velocity  field,  the  scalar  concentration  and  the  pressure.  The  dependence  of 
the  radial  velocity  on  the  actuation  parameters  is  pronounced,  while  the  variation  of  the  axial  and 
azimuthal  velocity  is  small  and  difficult  to  quantify.  We  have  therefore  maximized  the  integral  of 
the  radial  velocity  (in  the  cylindrical  system) 

f{x.,to)  =  fvV%c(r,e,(p,t0)dV'  (6) 

where  V  is  the  computational  domain  as  shown  in  Fig.  1.  The  integral  of  v^c  is  time  dependent.  A 
suitable  time  to  must  be  chosen  in  the  jet  simulation  for  the  evaluation  of  the  objective  function. 

Besides  (6)  we  have  tested  other  objective  functions.  The  distribution  of  the  passive  scalar 
corresponds  directly  to  the  amount  of  mixing  in  the  jet.  We  have  used  the  integral  of  the  scalar 
transported  to  the  outer  part  0  >  0q  of  the  computational  domain  as  an  objective  function 

/(x,i0)=  [  C(r,e,v,t0)  dV1  (7) 

where  C  is  the  concentration  of  the  scalar.  We  have  also  changed  (6)  by  limiting  the  integral  of 
vl  to  this  subdomain.  The  best  parameter  vector  produced  by  the  optimization  does  not  differ 
significantly  for  these  different  objective  functions.  On  the  other  hand,  if  we  calculate  the  integral 
of  v%  over  the  subdomain  10D  <  r  <  15 D,  the  result  differs  from  that  found  for  the  whole  domain. 
The  maximum  of  the  objective  function  value  is  now  located  at  lower  Strouhal  numbers.  For 
low  Strouhal  numbers  the  mode  that  determines  the  creation  of  vortices  saturates  and  is  replaced 
by  its  subharmonic  further  downstream.  The  first  pairing  of  vortices,  which  is  the  onset  of  the 
jet  spreading,  therefore  happens  further  away  from  the  orifice  for  small  Strouhal  numbers.  The 
corresponding  actuation  leads  to  a  lower  fitness  value  when  the  whole  domain  is  used  for  the 
integration  of  v% . 
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Another  possible  objective  function  is  the  centerline  velocity.  While  the  integral  of  v%  changes 
as  a  function  of  the  parameter  values  early  in  the  DNS,  the  centerline  velocity  does  not  show  clear 
variations  at  that  stage.  Different  objective  functions  have  been  investigated  by  Freund  and  Moin 
[8].  They  showed  that  volume  integrals  of  moments  of  the  scalar  concentration  and  integrals  of  the 
scalar  dissipation  are  suitable  metrics  to  quantify  mixing  in  a  jet. 

Our  computations  have  shown  that  the  optimization  strategy  tends  to  choose  at  least  one  of 
the  amplitudes  to  be  as  large  as  possible  within  the  given  limits.  Since  an  actuation  with  verj 
large  amplitudes  is  not  desirable,  we  have  kept  one  of  the  amplitudes,  usually  Aa,  constant  for  the 
dual-frequency  actuation  (3)  and  varied  the  ratio  of  amplitudes  Ah/Aa.  For  the  excitation  with 
one  frequency  (2),  we  have  introduced  a  penalty  function 

/(x)  =  [  vldV'-CAl  (8) 

which  seeks  to  maximize  the  spreading  of  the  jet  while  minimizing  the  amplitude  Ah.  The  constant 
C  >  0  is  chosen  such  that  the  two  terms  in  (8)  are  of  the  same  order  of  magnitude.  If  C  is  very 
large  the  second  term  dominates  the  objective  function  such  that  the  search  results  in  very  small 
amplitudes.  If  C  is  on  the  other  hand  very  small  the  effect  of  the  penalty  is  negligible. 

Shape  of  the  fitness  landscape 

The  performance  of  optimization  procedures  is  in  general  evaluated  by  applying  them  to  well  known 
test  functions.  In  order  to  provide  a  suitable  test  function,  we  have  chosen  the  case  where  the  two 
Strouhal  numbers  were  varied  while  the  amplitudes  were  kept  fixed.  We  were  able  to  numerically 
obtain  an  image  of  the  fitness  landscape  for  a  jet  at  low  Reynolds  number  Re  =  1500.  For  this 
purpose  we  calculated  the  fitness  value  on  an  equidistant  grid  in  (Sta,  /3)  space  (with  f3  =  Sta/ Sth). 
As  mentioned  earlier,  the  typical  CPU  time  for  a  full  jet  simulation  is  1200  node  hours.  Since  at 
least  200  evaluations  of  the  objective  function  are  necessary  to  identify  the  fitness  landscape,  the 
time  for  each  evaluation  needs  to  be  reduced  drastically. 

Figure  2  shows  how  the  objective  function  (6)  evolves  with  time  for  the  single- frequency  actuation 
(2).  The  three  curves  correspond  to  different  values  of  the  axial  Strouhal  number  at  constant 
amplitudes  Aa  =  0.025,  Ah  =  0.05  and  at  (3  =  2.  The  time  scale  is  the  normalized  scale  of  the  jet 
simulation.  (For  comparison,  the  time  scale  of  an  excitation  with  Sta  =  0.7  is  T  =  D/ (StV0)  =  0.23. 
The  time  for  the  DNS  of  a  fully  developed  jet  is  approximately  t  =  7.)  A  comparison  of  the  early 
stages  of  the  simulation  with  the  fully  developed  jets  has  shown  that  pronounced  jet  spreading 
corresponds  to  a  large  value  of  the  objective  function  at  these  early  times.  We  have  therefore 
determined  the  fitness  value  of  the  parameter  vector  at  time  to  =  1.6  where  the  objective  function 
values  already  clearly  differ.  Fig.  2  shows  the  objective  function  for  the  Strouhal  number  Sta  =  0.7 
that  leads  to  a  high  value  at  t  =  1.6  and  for  Sta  =  0.2  and  Sta  =  0.9  that  result  in  much  lower  values. 
The  time  scale  shown  in  Fig.  2  is  not  proportional  to  the  CPU  time  needed  for  the  simulation. 
Note  that  the  time  integration  of  the  DNS  proceeds  much  faster  at  the  beginning  of  the  simulation, 
starting  from  a  laminar  flow.  Calculating  only  the  early  stages  of  the  jet  development  therefore 
leads  to  a  large  increase  of  the  optimization  speed. 

In  order  to  decrease  the  computational  time,  we  have  used  a  coarse  grid  for  the  optimization 
procedure.  A  comparison  of  calculations  on  different  grids  has  shown  that  most  of  the  small 
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scale  structures  are  resolved  on  a  128  x  96  x  64  or  finer  grid.  The  overall  behavior  of  the  jet,  in 
particular  at  early  stages  of  the  simulation,  however  is  captured  by  simulations  on  a  coarser  grid. 
The  objective  function  shows  a  very  similar  behavior  for  early  time  t  <to  on  coarse  and  fine  grids. 
We  note  that  under-resolved  simulations  can  only  be  used  for  a  rough  estimate  of  the  jet  dynamics, 
but  since  the  jet  flow  is  controlled  by  large  scale  structures,  the  overall  dynamic  is  preserved  b\ 
this  approximation.  Fig.  3  shows  sketches  of  the  fitness  landscape,  which  have  been  obtained  by 
calculating  the  objective  function  value,  i.e.  the  value  of  (6)  at  time  t  -  t0,  for  an  equidistant,  grid 
in  {Sta,{3) —space.  This  calculation  has  been  performed  in  order  to  compare  the  fitness  landscape 
for  different  grid  size  and  to  provide  a  test  case  for  the  optimization  schemes. 

In  Fig.  3  a),  the  DNS  was  performed  on  a  128  x  96  x  64  grid.  The  landscape  shows  a 
pronounced  global  maximum,  which  appears  at  the  parameter  values  {Sta,0)  «  (0.68,2.0),  as  well 
as  local  maxima  with  lower  amplitude  at  (Sta,(3)  ~  (0.9, 2.8)  and  ( Sta,(3 )  ~  (1.1, 3. o).  These 
maxima  are  located  on  a  straight  line  /?  »  3  •  St0,  which  corresponds  to  a  helical  Strouhal  number 
Sth  «  0.33.  In  fact,  most  fitness  values  along  this  line  are  significantly  higher  than  those  in  their 
vicinity.  Another  range  of  higher  values  can  be  seen  for  Sta  a  0.5  and  3  2.  For  Fig.  3  b), 

a  coarse  grid  64  x  64  x  32  has  been  chosen.  Although  the  extrema  are  not  as  pronounced  as  in 
figure  a),  the  landscape  roughly  has  the  same  shape  for  simulation  on  this  coarse  grid.  The  global 
optimum  appears  again  at  ( Sta,/3 )  a  (0.7, 2.0)  and  there  is  a  range  of  high  values  along  the  line 
Sth  a  0.3.  Using  the  approximations  described  in  this  section,  the  CPU  time  necessary  for  one 
evaluation  of  the  objective  function  is  approximately  3  node  hours  on  a  SGI-ORIGIN  2000.  In  fact, 
the  optimization  leads  to  good  results  while  keeping  the  computational  time  within  reasonable 
limits. 

In  the  next  section,  we  will  use  our  knowledge  of  the  fitness  landscape  to  compare  the  per¬ 
formance  of  the  different  optimization  procedures.  Since  the  numerical  estimation  of  the  fitness 
landscape  is  computationally  expensive,  we  have  not  repeated  this  procedure  for  other  cases  that 
were  optimized  -  different  parameters,  different  objective  functions  and  higher  Reynolds  numbers. 
Instead,  we  assume  that  the  test  case  chosen  is  representative  for  our  optimization  procedure. 


Results  of  the  optimization 

In  this  section  we  describe  the  results  obtained  by  optimizing  the  objective  function  (6)  for  different 
actuation  parameters.  We  present  results  for  actuation  using  helical  actuation  at  the  jet  orifice. 
Since  we  only  varied  one  parameter,  we  used  the  simple  (1  +  1)  evolution  strategy.  A  combined 
axial  and  helical  forcing  has  been  optimized  with  simulated  annealing  and  with  a  multi-membered 
evolution  strategy.  Using  our  knowledge  of  the  fitness  landscape  described  in  the  previous  section, 
the  performance  of  the  optimization  strategies  is  compared. 

One-frequency  actuation 

First,  we  varied  only  the  Strouhal  number  of  the  actuation  (2)  and  kept  the  amplitude  fixed  at 
Ah  =  0.05.  The  optimization  was  done  with  the  objective  function  (6)  and  the  (1  +  1)  evolution 
strategy.  Starting  from  an  initial  value  Sth  =  0-5,  the  evolution  path  was  found  to  approach  an 
optimum  in  the  27th  generation  (28  evaluations  of  the  fitness  function)  at  St  =  0.36.  We  repeated 
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the  calculation  with  an  objective  function  that  integrates  v~  only  in  the  subdomain  10D  <r  <  15  D. 
This  neglects  spreading  in  the  part  of  the  computational  domain  close  to  the  orifice  and.  therefore, 
favors  lower  Strouhal  numbers.  As  a  result  of  the  optimization,  we  found  St  =  0.28.  While  Parekh 
et  al  found  the  most  intense  motion  of  the  jet  for  St  ~  0.2,  they  also  observed  efficient  mixing  at 
St  »  0.27,  which  is  close  to  our  result,  and  at  St  ^  0.4. 

As  mentioned  in  the  previous  section,  the  simultaneous  optimization  of  Strouhal  number  and 
amplitude  Ah  results  in  the  maximum  possible  amplitude  within  the  given  limits.  In  order  to 
keep  the  amplitude  low  we  repeated  the  optimization  with  the  function  (8)  that  penalizes  large 
amplitudes.  The  factor  describing  the  weight  of  the  penalty  was  chosen  as  C  =  5  and  the  allowed 
range  for  the  amplitude  was  0.02  <  Ah  5;  0.12.  Starting  with  the  initial  parameter  vector  x 
(0.37,0.05)  for  Strouhal  number  and  amplitude,  we  found  the  maximum  fitness  value  at  x  = 
(0.29,0.08).  If  we  choose  a  larger  factor  C,  smaller  amplitudes  are  preferred.  For  C  =  15,  the 
optimum  is  found  at  x  =  (0.30,0.045).  The  Strouhal  number  is  approximately  the  same  for  both 
values  of  C.  In  the  literature  of  evolution  programs  sophisticated  ways  of  dealing  with  multi¬ 
objective  optimization  problems  have  been  proposed  [35].  However,  we  have  chosen  the  objective 
function  as  shown  in  Eqn.  (8)  for  simplicity. 

For  the  results  described  so  far,  the  jet  simulation  was  started  at  time  t  =  0  from  the  initial 
laminar  flow  at  the  orifice.  The  optimization  has  been  repeated  with  an  initial  non-laminar  velocity 
field,  which  was  obtained  from  a  jet  simulation  with  small  axial  and  random  forcing.  The  parameter 
vector  found  by  the  evolution  strategy  did  not  differ  significantly. 

Two-frequency  actuation 

It  has  been  shown  in  many  experiments  that  the  combination  of  axial  and  helical  actuation  can 
lead  to  a  high  mixture  fraction  in  the  jet  [2,  3,  9].  We  will  now  present  results  that  we  obtained  by 
optimizing  the  parameters  of  the  dual- frequency  actuation  (3).  Although  this  actuation  contains 
four  parameters,  we  varied  only  three  of  them  in  order  to  reduce  the  dimension  of  the  search  space. 
For  most  runs  we  kept  the  axial  amplitude  constant  at  Aa  =  0.025  and  optimized  the  parameter 
vector  x  =  (Sta,  (3,Ah),  /3  =  Sta/Sth.  Since  we  expected  the  search  to  be  more  difficult  than  for  the 
variation  of  just  one  parameter,  we  used  a  (/x  +  A)  strategy  with  fi  =  2  parents  and  A  =  5  offspring 
and  an  adaption  of  the  step  length  with  the  covariance  matrix  adaption  scheme.  Although  this 
still  corresponds  to  a  small  population,  it  increases  the  efficiency  of  the  search  while  keeping  the 
CPU  time  needed  for  the  evaluation  of  each  generation  within  reasonable  limits.  The  limits  for 
the  parameters  were  chosen  as  described  in  section  2.  In  particular,  we  have  set  the  upper  limit 
for  the  helical  amplitude  to  Ah  =  0.075.  Starting  from  the  parameter  vector  x0,i  =  (0.7,2.2,0.05) 
which  corresponds  to  an  objective  function  value  /(x o,i)  =  2.19  the  best  parameter  vector  X(,e5<  = 
(0.66, 2.1,0.075)  with  fitness  value  finest)  =  3.12  was  obtained  after  25  generations.  We  repeated 
the  simulation  with  different  limits  for  the  helical  amplitude  0.0075  <  Ah  <  0.075  and  a  different 
start  vector  x0,2  =  (0.5,2.5,0.0175)  with  /(x0,2)  =  1-38,  and  reached  the  same  optimum  after  40 
generation  which  corresponds  to  201  evaluations  of  the  objective  function.  For  comparison,  the 
objective  function  value  of  a  jet  simulation  without  actuation  (laminar  flow)  is  /(0)  =  0.50). 

We  have  calculated  the  objective  function  value  in  the  vicinity  of  the  global  optimum  to  estimate 
the  sensitivity  of  the  result  to  changes  of  the  parameter  vectors.  We  found  that  the  objective 
function  value  drops  to  98%  if  the  parameters  are  changed  by  A Sta  =  0.01  and  A/3  =  0.15  (The 
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value  of  98%  has  been  calculated  with  respect  to  the  total  range  of  the  objective  function  value 
found  in  our  optimization,  /(x)  €  [0.50,3.12].)  Therefore  errors  larger  than  a  few  percent  should 
be  avoided  in  applications. 

The  results  of  the  optimization  described  so  far  suggest  that  the  search  strategy  always  chooses 
the  maximal  allowed  helical  amplitude.  While  the  objective  function  value  grows  with  the  ampli¬ 
tude,  variations  of  Ah  do  not  have  a  large  influence  on  the  shape  of  the  fitness  landscape.  In  order 
to  confirm  this  assumption  we  have  repeated  the  calculation  with  the  parameter  vector  ( Sta,@ ) 
and  kept  the  amplitudes  fixed  at  Aa  =  0.025,  Ah, i  =  0.075  and  Ah, 2  =  0-05  respectively.  We  found 
the  global  optima  at  x6estii  =  (0.65,2.0)  and  x6es<i2  =  (0.66,2.1),  in  agreement  with  the  previous 
results.  This  confirms  that  the  frequencies  are  the  relevant  parameters  for  the  maximization  of  jet 
spreading  for  the  amplitude  range  considered  in  this  study.  However,  it  has  been  shown  in  [28]  that 
for  axial  forcing  Strouhal  numbers  that  enhance  mixing  are  different  at  very  low  {A  <  0.01)  and  at 
high  forcing  levels. 

In  order  to  test  the  interdependence  of  the  two  actuation  frequencies,  we  tried  to  optimize  the 
objective  function  separately  for  the  two  frequencies.  We  found  that  the  location  of  the  maximum 
differs  in  this  case  and  that  the  corresponding  fitness  value  is  lower.  The  Strouhal  numbers  are  thus 
not  independent  and  should  be  optimized  simultaneously.  This  result  is  not  surprising  because  the 
spreading  pattern  of  the  jet  is  due  to  the  interaction  of  the  various  modes  in  the  jet. 

Apart  from  the  function  (6)  we  used  the  integral  of  the  scalar  concentration  (7)  as  objective 
function.  We  also  changed  the  objective  function  by  extending  the  integral  of  v%  or  the  scalar  not 
to  the  whole  domain,  but  only  to  the  outer  part  6  >  Qq.  For  all  these  choices  the  same  optimal 
Strouhal  numbers  ( Sta,/3 )  «  (0.66,2.1)  were  found.  This  shows  that  our  optimization  does  not 
depend  critically  on  the  choice  of  the  quantity  (radial  velocity,  scalar  concentration)  measured  by 
the  objective  function. 

Comparison  of  the  optimization  procedures 

In  order  to  compare  the  different  search  strategies  for  our  test  case,  the  optimization  of  the  objective 
function  (6)  for  the  dual- frequency  actuation  (3)  with  two  parameters  xo  =  (Sta,f3)  and  fixed 
amplitudes,  we  have  started  the  optimization  at  different  points  in  phase  space.  We  determined 
whether  the  global  optimum  was  found  and  the  number  of  fitness  function  evaluations  necessary  to 
reach  it.  The  optimization  was  done  for  our  test  case,  the  dual- frequency  actuation  (3)  with  two 
parameters  xo  =  ( Sta,j3 ).  The  amplitudes  were  kept  constant  at  Aa  =  0.025  and  A/,  =  0.05  and 
the  fitness  value  was  calculated  at  to  =  1-6. 

The  starting  points  for  the  optimization  procedure  were  chosen  from  a  contour  plot  of  the 
fitness  landscape,  shown  in  Fig.  4.  They  were  close  to  the  global  maximum,  close  to  a  local 
maximum  and  far  away  from  both.  Using  simulated  annealing,  the  (1  +  1)  and  the  (2,5)  evolution 
strategies,  we  have  recorded  the  number  of  objective  function  evaluations  necessary  to  reach  the 
global  optimum.  This  number  determines  the  required  CPU  time  which  is  the  limiting  factor  for 
our  calculations.  We  assume  that  the  optimization  has  converged  when  the  final  fitness  value  is  98% 
of  the  maximal  value  at  the  global  optimum  or  higher  and  if  the  actuation  parameters  are  within 
the  corresponding  range  as  described  in  Sec.  .  The  results  of  the  optimization  are  summarized  in 
table  1.  In  all  cases,  the  global  optimum  was  found.  The  listed  numbers  are  lower  than  the  total 
number  of  function  evaluations  carried  out  for  each  run:  When  the  optimum  has  been  reached  the 
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procedure  continues  searching  for  better  solutions  and  stops  only  if  no  improvement  has  been  made 
for  a  large  number  of  iterations.  Note  that  for  the  (2, 5)  strategy  each  generation  corresponds  to 
5  evaluations  of  the  objective  function.  On  average,  simulated  annealing  and  the  (2,5)  evolution 
strategy  with  covariance  matrix  adaption  of  the  steplength  were  able  to  find  the  global  optimum 
using  a  comparable  number  of  function  evaluations,  the  number  being  slightly  lower  for  simulated 
annealing.  The  (1  +  1)  strategy  on  the  other  hand  needed  in  general  a  larger  number  of  steps.  The 
fact  that  the  number  of  function  evaluations  for  the  (2, 5)  evolution  strategy'  is  not  lower  than  that 
for  simulated  annealing  indicates  that  the  number  of  search  steps  that  can  be  achieved  for  this 
computationally  expensive  problem  is  too  small  to  allow  efficient  use  of  the  information  gathered 
on  the  search  path.  We  have  therefore  used  this  method  mainly  for  the  less  expensive  optimization 
at  low  Reynolds  numbers  and  simulated  annealing  for  the  optimization  at  higher  Re  which  will  be 
presented  in  the  next  section. 

For  the  start  location  xo  =  (0.65, 1.6)  close  to  the  global  optimum,  the  (2,5)  ES  and  simulated 
annealing  required  56  and  60  objective  function  evaluations.  The  reason  for  this  is  that  a  large  initial 
step  size  was  chosen  for  all  strategies  such  that  the  vicinity  of  the  initial  point  is  not  necessarily 
searched  first.  For  the  start  location  xo  =  (0.5, 2.5),  both  evolution  strategies  had  difficulties  finding 
the  global  optimum,  requiring  201  fitness  function  evaluations  even  for  the  (2,5)  ES.  Simulated 
annealing  on  the  other  hand  needed  102  function  evaluations  when  it  was  started  from  xo  = 
(0.9, 2.8),  in  contrast  to  31  evaluations  for  the  (2,5)  ES. 

Figure  5  gives  an  example  for  the  convergence  of  the  three  optimization  strategies.  For  the  (2, 5) 
ES,  the  fitness  value  of  the  best  offspring  is  shown  as  a  function  of  the  generation  number.  The 
resulting  curve  is  not  monotonic  because  the  offspring  are  obtained  by  random  mutation  of  the 
average  parent  vectors  and  therefore  may  have  lower  fitness  values.  For  simulated  annealing  and 
the  (1  +  1)  ES  on  the  other  hand,  the  best  solution  is  kept  until  it  is  replaced  by  a  better  one.  the 
best  value  found  so  far  is  shown  as  a  function  of  the  generation  number. 


Xo 

(0.45,3.1) 

(0.5, 2.5) 

(0.65, 1.6) 

(0.9, 2.8) 

(1.0, 2.0) 

(1.0, 3.2) 

sim.  annealing 

52 

42 

60 

102 

41 

25 

(2,5)  ES 

61 

201 

56 

31 

41 

51  • 

(1+1)  ES 

153 

262 

58 

111 

71 

62 

Table  1:  Comparison  of  the  optimization  strategies:  number  of  objective  function  evaluations 
required  to  reach  the  global  optimum 

Simulation  at  higher  Reynolds  number  using  LES 

For  many  applications  of  turbulent  jet  control,  the  Reynolds  number  of  the  flow  is  much  higher  than 
that  used  for  the  simulation  presented  so  far.  For  the  control  of  jets  at  higher  Reynolds  numbers 
we  included  a  subgrid  scale  model  in  our  simulation  to  keep  the  computational  cost  affordable.  The 
idea  of  large  eddy  simulation  (LES)  is  to  resolve  only  the  large  scales  of  the  flow  while  modeling  the 
contribution  of  the  small  (subgrid)  scales.  The  LES  equations  are  obtained  by  applying  a  spatial 
filter  to  the  momentum  equations.  In  order  to  account  for  the  subgrid-scale  contribution  to  the  flux 
of  momentum  the  molecular  viscosity  is  augmented  by  an  eddy  viscosity  ut  =  CA2|S|,  where  S  is 
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the  subgrid  strain  tensor  and  C  the  Smagor insky  constant.  We  have  used  a  large  eddy  simulation 
(LES)  based  on  the  dynamic  procedure  by  Germano  at  al.  [36],  which  determines  the  constant  C 
as  a  function  of  space  and  time.  The  dynamic  procedure  does  not  contain  any  model  constants. 

We  have  used  the  combination  of  optimization  algorithm  and  LES  to  search  for  actuation  that 
maximizes  mixing  in  a  Re  =  6000  turbulent  jet.  The  simulations  are  more  expensive  than  the  DNS 
at  Re  —  1500  such  that  a  numerical  estimate  of  the  fitness  landscape  is  impossible.  Instead,  we  use 
the  conclusions  drawn  from  our  test  case:  We  confine  the  search  to  the  space  of  the  two  Strouhal 
numbers,  assuming  that  they  are  the  relevant  parameters  for  control  of  jet  mixing.  The  parameters 
were  varied  in  the  range  Sta  G  [0.2, 1.2]  and  /?  G  [1.6, 3.2]  and  the  amplitudes  were  kept  constant 
at  Aa  =  0.025  and  Ah  =  0.075.  The  optimization  was  performed  with  the  simulated  annealing 
algorithm,  which  on  average  needed  the  lowest  number  of  objective  function  evaluations  for  the 
test  case. 

First,  we  compared  simulations  on  grids  with  different  resolution.  We  found  that  a  grid  with 
160  x  120  x  32  points  in  the  radial,  tangential  and  azimuthal  direction  respectively  is  small  enough  to 
allow  repeated  calculation  of  the  objective  function,  while  being  sufficiently  large  to  give  a  reliable 
estimate  of  the  jet  dynamics.  Comparing  the  objective  function  (6)  as  a  function  of  time  for  different 
parameter  values,  we  have  chosen  the  time  for  the  evaluation  of  the  function  as  to  =  1-8.  The  CPU 
time  for  one  evaluation  of  the  objective  function  was  32  node  hours,  as  compared  to  1800  node 
hours  for  the  full  LES  on  a  192  x  128  x  64  grid. 

Starting  the  search  at  xo  =  ( Sta,(3 )  =  (0.66, 2.1)  (the  optimum  for  Re  —  1500),  the  best  value 
was  reached  at  x  =  (0.79,  2.18)  which  corresponds  to  ( Sta ,  Sth)  =  (0.79,0.36)  after  52  evaluations 
of  the  objective  function.  Comparing  the  absolute  fitness  values  we  found  that  the  maximum 
value  reached  at  this  Reynolds  number  is  approximately  1.5  times  higher  than  that  for  Re  =  1500, 
suggesting  that  the  spreading  of  the  jet  is  larger  for  the  Re  =  6000  jet.  Apart  from  the  maximum 
at  x  =  (0.79,2.18)  with  fitness  value  /  =  5.0  we  found  parameters  with  comparably  high  value  at 
x  =  (0.90, 2.0)  with  /(x)  =  4.9  and  x  =  (0.58, 2.6)  with  /(x)  =  4.8.  The  lowest  value  found  was 
/  =  3.8.  In  general  the  objective  function  values  were  higher  for  Sta  >  0.5  than  for  lower  Strouhal 
numbers.  Apparently  the  fitness  landscape  has  a  structure  that  is  different  from  that  found  for 
Re  =  1500.  In  particular,  high  fitness  values  do  not  only  appear  at  Sth  ~  0-33.  Since  we  are 
only  able  to  determine  a  small  number  of  points  on  the  fitness  landscape,  a  guess  of  its  structure 
is  difficult.  Like  in  most  optimization  problems  there  is  no  certainty  about  whether  the  global 
optimum  has  been  reached.  However,  we  have  found  a  good  solution  within  reasonable  time. 


Discussion  of  the  results 

For  the  best  actuation  parameters  found  by  the  evolution  strategies  we  have  repeated  the  DNS 
of  the  jet  on  the  fine  grid  described  in  the  second  section.  Figs.  6  and  7  show  the  passive  scalar 
concentration  obtained  when  the  helical  actuation  (2)  is  applied  at  the  orifice.  For  the  best  Strouhal 
number  Sth  ~  0-36,  Fig.  6  shows  snapshots  taken  at  time  t  =  9  (on  the  normalized  time  scale  of 
the  jet).  Different  shades  of  grey  denote  different  concentration  C  of  the  scalar.  The  concentration 
is  approximately  one  in  the  inner  (dark)  region  and  zero  far  outside  (white  region).  The  figure 
shows  the  jet  in  the  plane  of  the  actuation,  <p  =  0  (left),  and  in  the  plane  tp  —  7t/2  (right).  The  jet 
spreads  rapidly  in  the  plane  of  the  actuation  and  contracts  in  the  orthogonal  plane.  Although  the 
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amplitude  of  the  actuation  is  small  (Ah  =  0.05),  the  jet  spreads  at  a  large  angle  and  the  jet  column 
shows  a  strong  flapping  motion.  The  jet  column  disintegrates  towards  the  end  of  the  computational 
domain,  but  regions  of  high  concentration  C  w  1  remain  near  the  centerline  of  the  jet.  Fig.  7  shows 
the  best  result  obtained  by  maximizing  the  spreading  in  half  of  the  computational  domain.  Since 
the  optimal  Strouhal  number  is  smaller  in  this  case,  the  jet  spreads  further  downstream  from  the 
orifice  than  in  the  previous  case. 

Figure  8  shows  a  snapshot  of  the  jet  actuated  with  the  parameters  Xbest  =  {Sta,Sth.  Ah)  = 
(0.66,0.31,0.075)  and  Aa  =  0.025.  This  dual-frequency  actuation  leads  to  an  impressive  spreading 
of  the  jet.  The  center  of  the  jet  shows  a  strong  flapping  motion,  which  is  due  to  the  large  amplitude 
of  the  helical  forcing.  Comparing  Figs.  8,  6  and  7  we  find  that  the  shape  of  the  jet  differs  for  the 
three  types  of  actuations.  The  overall  spreading,  i.e.  the  area,  to  which  a  significant  amount  of 
the  tracer  is  transported,  is  comparable,  although  the  total  amplitude  is  smaller  for  the  jets  with 
single-frequency  forcing.  The  flapping  motion  of  the  jet  column  however  is  most  pronounced  for 
the  dual-frequency  actuation.  The  jet  column  is  completely  dispersed  near  the  end  of  the  domain, 
indicating  that  the  jet  bifurcates.  Although  a  clear  spreading  in  two  branches  is  not  visible,  there 
are  areas  on  both  sides  of  the  centerline  with  much  higher  concentration  of  the  scalar  and  much 
lower  radial  velocity  (not  shown)  than  on  the  centerline.  For  comparison  we  also  show  a  natural 
jet,  calculated  for  a  slightly  higher  Reynolds  number,  in  Fig.  9.  A  larger  domain  is  chosen  since 
the  potential  core  is  much  longer  in  this  case.  The  spreading  angle  of  the  jet  is  approximately  10°. 
A  quantitative  comparison  of  the  jets  obtained  by  our  simulation  with  the  experimental  results  by 
Parekh  et  al.  (1996)  is  difficult  because  the  Reynolds  numbers  differ.  In  addition,  different  types 
of  perturbation  are  used  in  the  experiment  and  simulation,  and  hence  spreading  angles  will  always 
be  slightly  different. 

In  Fig.  10  a  snapshot  of  the  scalar  concentration  is  plotted  for  a  jet  at  Re  —  6000  which  was 
obtained  using  the  actuation  with  the  best  Strouhal  numbers  Sta  =  0.79  and  Sth  =  0.36  found  for 
this  Reynolds  number.  The  flow  was  calculated  on  a  grid  with  192  x  160  x  64  points.  Again,  a 
strong  flapping  motion  of  the  jet  column  is  visible.  Apart  from  the  large  spreading  of  the  jet  it  can 
be  seen  that  the  jet  column  disintegrates  soon  after  the  initial  pairing.  Further  downstream  the 
values  of  the  scalar  concentration  are  low,  indicating  good  mixing  in  the  jet. 

The  streamwise  velocity  and  concentration  of  the  scalar  on  the  centerline  of  the  jet  is  evaluated 
in  order  to  obtain  a  more  quantitative  measure  of  turbulent  mixing.  It  is  known  that  the  decay 
rate  of  round  jets  is  proportional  to  [z  -  z0]~l  where  zq  is  the  virtual  origin  of  the  jet  and  z  the 
distance  measured  from  the  jet  orifice  [37].  The  perturbations  given  by  equations  (2)  and  (3)  axe 
not  axisymmetric  and  therefore  linear  decay  of  the  jet  is  not  expected.  In  Figure  11  (left)  we  show 
the  centerline  velocity  obtained  from  the  flapping  and  bifurcating  jet  computations.  For  comparison 
we  include  the  profile  for  a  round  axisymmetric  jet  (standard  jet)  at  comparable  Reynolds  number. 
For  both  the  flapping  and  bifurcating  jets  the  centerline  velocity  starts  to  drop  near  z/D  =  5,  which 
is  much  earlier  than  in  the  standard  jet.  The  decay  rate  (slope  of  the  curve  in  Fig.  11  (right)) 
of  the  flapping  jet  is  comparable  to  that  of  the  standard  jet  while  the  bifurcating  jet  decays  much 
faster.  For  the  flapping  jet  the  decay  is  approximately  linear,  for  the  bifurcating  jet  superlinear. 

In  our  simulation  it  can  be  assumed  that  for  values  of  the  scalar  concentration  C  ~  0.5  the 
flow  is  well  mixed.  In  Fig.  12  we  show  the  centerline  scalar  concentration  for  the  bifurcating,  the 
flapping  and  the  natural  jets.  The  actuation  leads  to  a  much  earlier  decay  of  concentration  than 
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for  the  natural  jet.  For  the  flapping  actuation,  the  concentration  reaches  C  a  0.4  on  the  centerline 
at  the  end  of  the  domain.  For  the  bifurcating  jet,  the  decay  starts  even  earlier  and  the  mixing  is 
very  efficient.  Towards  the  end  of  the  domain,  the  average  concentration  reaches  low  values  due  to 
the  bifurcation  of  the  jet,  which  transports  the  scalar  away  from  the  centerline. 

Fig.  13  shows  the  total  turbulent  flux  of  the  scalar  \Ju'TC'  +u'zC'  in  the  plane  of  the  actuation 
for  both  the  flapping  and  bifurcating  jet.  For  the  flapping  jet  the  flux  has  high  values  near  the 
centerline.  The  bifurcating  actuation  on  the  other  hand  directs  the  total  flux  outwards.  The 
bifurcating  jet  therefore  has  by  far  the  best  mixing  properties.  From  the  figures  presented  in  this 
section  it  is  clear  that  there  is  a  certain  amount  of  scatter  in  the  DNS  data  which  is  due  to  the 
limited  size  of  the  statistical  sample. 

The  main  results  of  our  optimization  are  summarized  in  table  2.  In  all  cases,  the  spreading  of 
the  jet  is  most  pronounced  at  helical  Strouhal  numbers  St  ~  0.3. ..0.36.  It  differs  from  the  natural 
Strouhal  number,  for  which  linear  instability  analysis  predicts  the  largest  amplification  of  signals 
(Michalke  1984),  but  is  in  agreement  with  the  preferred  frequency  of  the  jet  found  by  Crow  et  al. 
[26]  and  Mankbadi  [28].  The  optimal  axial  Strouhal  number  found  by  the  optimization  procedure 
is  approximately  twice  as  large  as  the  preferred  frequency,  as  predicted  by  Mankbadi  [28].  The 
Strouhal  numbers  found  are  in  very  good  agreement  with  experimental  results  by  Lee  at  al.  [2] . 

In  the  optimization  procedure  described  in  this  paper  some  arbitrary  choices  have  been  made 
that  may  influence  the  results.  The  optimal  parameter  vector  depends  on  the  choice  of  the  time  to, 
at  which  the  objective  function  value  is  determined.  It  also  depends  on  the  domain  that  is  chosen 
for  the  integration  of  vj  or  in  general  on  the  length  of  the  computational  domain.  However,  for 
each  choice  of  the  domain  actuation  parameters  were  found  that  lead  to  a  large  spreading  of  the 
jet. 


Re  =  1500 


Sta 

sth 

Aa 

Aa 

single-freq.  actuation,  eq.  (2) 

0.36 

- 

0.05  (fixed) 

single-freq.  actuation,  penalty  C  —  5 

0.29 

- 

0.08 

two-freq.  actuation,  eq.  (3) 

0.66 

0.31 

0.025  (fixed) 

0.075 

two-freq.  actuation,  separate  opt.  of  Sta,Sth 

0.72 

0.29 

0.025  (fixed) 

0.075  (fixed) 

Re  =  6000 


Sta 

sth 

Aa 

Aa 

two-freq.  actuation,  eq.  (3) 

0.79 

0.36 

0.025  (fixed) 

0.075  (fixed) 

Table  2:  Best  actuation  parameters  found  by  the  optimization  strategies. 


Conclusions 

In  this  paper  we  have  combined  stochastic  search  programs  and  numerical  simulation  of  a  jet  to 
search  for  actuation  parameters  (Strouhal  numbers  and  amplitudes)  that  enhance  mixing  in  the 
jet.  We  have  used  different  objective  functions  that  evaluate  the  performance  of  a  given  actuation 
and  found  that  the  integral  of  the  squared  radial  velocity  and  the  integral  of  the  concentration 
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of  a  passive  scalar  are  suitable  measures.  Since  jet  flows  are  governed  by  the  dynamics  of  the 
large  eddies,  we  have  used  a  coarse  grid  for  the  DNS  and  LES  in  the  optimization  procedure.  We 
have  found  that  the  spreading  of  the  jet  and  hence  the  amount  of  mixing  depends  very  much  on 
the  actuation  Strouhal  numbers  and  that  it  grows  with  the  amplitudes.  The  main  result  is  that  a 
bifurcating  perturbation  gives  better  mixing  and  a  faster  decay  of  the  center  line  velocity  and  scalar 
concentration  than  a  flapping  perturbation.  We  found  the  most  pronounced  spreading  for  helical 
and  axial  Strouhal  numbers  Sth  «  0.3... 0.36,  Sta  «  2  •  Sth  at  Reynolds  number  Re  =  1500  and 
Re  =  6000.  The  spreading  of  the  jet  as  measured  by  our  objective  function  is  larger  at  the  higher 
Reynolds  number.  For  applications,  it  is  necessary  to  have  a  large  spreading  not  only  in  one  plane 
but  in  the  whole  three  dimensional  domain.  This  can  be  obtained  by  replacing  the  phase- locked 
helical  forcing  used  in  our  simulation  by  a  helical  actuation  that  is  rotating  around  the  orifice. 

Stochastic  optimization  techniques  appear  to  be  very  efficient  for  optimization  of  the  jet  ac¬ 
tuation.  Comparing  different  methods  we  have  found  that  multi-membered  evolution  strategies 
and  in  particular  simulated  annealing  are  well  suited  for  this  problem.  Although  convergence  to 
a  local  optimum  can  never  be  completely  excluded,  it  can  be  stated  that  our  strategy  is  able  to 
find  a  good  solution  within  reasonable  time.  In  a  real  experiment  the  perturbation  is  not  described 
by  a  simple  mathematical  function  like  in  our  simulation,  and  therefore  a  one  to  one  comparison 
between  simulation  and  experiment  is  difficult.  However,  we  expect  that  the  frequency  range  found 
to  be  best  in  our  simulations  is  similar  to  that  for  experiments  under  comparable  conditions.  We 
also  expect  that  stochastic  search  strategies  can  be  used  to  optimize  physical  experiments.  Future 
work  will  concentrate  on  the  use  of  more  realistic  actuators  and  on  simulations  at  higher  Reynolds 
numbers. 
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Figure  2:  Objective  function  /  vs.  time  for  a  jet  forced  with  dual¬ 
frequency  excitation  (2),  Aa  =  0.025,  Ah  =  0.05  and  /3  =  2; 
the  fitness  value  is  determined  at  t  =  1.6  (128  x  96  x  64  grid) 


Figure  3:  Fitness  landscape  for  simulation  on  a  a)  128  x  96  x  64  and  b)  64  x  64  x  32  grid  (fitness 
value  in  arbitraty  units),  Aa  =  0.025,  Ah  =  0.05. 


Sta 

Figure  4:  Contour  plot  of  the  fitness  landscape  for  simulation  on  the  128  x  96  x  64  grid,  same 
parameter  values  as  in  Fig.  3 
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Figure  5:  Convergence  of  the  evolution  strategies  and  simulated  annealing  (SA)  for  the  parameter 
vector  x  =  ( Sta,/3 ),  xq  =  (0.45, 3.1),  and  Aa  =  0.025,  Ah  =  0.075  fixed. 
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Figure  6:  Snapshots  of  the  passive  scalar  concentration  at  time  t  =  9  for  a  jet  excited  with  the 
single-frequency  actuation  (2)  in  the  plane  of  the  actuation  ip  =  0  (left)  and  the  plane  tp  =  -n/2 
(right);  St  =  0.036  and  Ah  =  0.05;  Re  =  1500. 


Figure  7:  Snapshot  of  the  scalar  concentration  at  time  t  =  7  for  a  jet  actuated  with  single-frequency 
actuation  (2),  St  =  0.028  and  Ah  =  0.05;  Re  =  BSQ0. 


Figure  8:  Snapshot  of  the  scalar  concentration  at  time  t  —  7  for  a  jet  actuated  with  dual-frequency 
forcing  (3)  and  Sta  =  0.66,  Sth  =  0.31,  Aa  =  0.025,  Ah  =  0.075;  Re  =  1500. 


Figure  9:  Snapshot  of  the  scalar  concentration  of  a  natural  jet  without  significant  actuation,  Re  = 
2000. 


Figure  10:  Snapshot  of  the  scalar  concentration  at  time  t  =  8  for  a  jet  at  Re  =  6000,  Sta  - 
0.79,  Sth  =  0.36,  Aa  =  0.025  and  Ah  =  0.075. 


Figure  11:  Left:  The  averaged  centerline  velocity  for  the  bifurcating,  flapping  and  standard  jet 
Right:  The  inverse  of  the  centerline  velocity  for  the  flapping,  bifurcating  and  standard  jet. 


Appendix  D 


Instability  Modes  of  Cylindrical  Jet 
Based  on  Incompressible  Inviscid  Model: 
Piecewise  Continuous  Velocity  Profile 


A.  Pal 
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Instability  Modes  of  Cylindrical  Jet 
Based  on  Incompressible,  Inviscid  Model; 
Piecewise  Continuous  Velocity  Profile 

Alex  Pal1 


Stability  analysis  of  shear  layers  of  parallel,  incompressible,  inviscid  flows  is 
discussed  in  Drazin  and  Reid.2  For  the  solution  of  the  resulting  eigenvalue  problem  a 
numerical  method  is  proposed  by  assuming  a  piecewise  linear  velocity  profile.  Even 
though  the  analysis  assumes  temporal  instability,  the  method  can  be  applied  to  the  spatial 
instability  problem  as  well.  This  memorandum  proposes  an  analogous  method  for  the 
axisymmetric  jet. 

In  the  axisymmetric  case  the  stream  function  in  the  individual  subregions  of  the 
flow  (in  the  jet  potential  core,  in  the  shear  layer,  outside  the  shear  layer)  can  be  expressed 

in  closed  form  just  as  in  the  parallel  flow  case. 

The  resulting  closed  form  dispersion  equation  connecting  the  normalized 
frequency  u  and  wave  number  a  is  in  terms  of  modified  Bessel  functions  of  the  radius. 

The  use  of  these  functions  and  the  numerical  approach  recommended  are  discussed  in  the 
second  half  of  this  report. 

Governing  equations 

Similar  to  the  planar  case,  discussed  in  Drazin  and  Reid,  the  following  assumptions 
are  made: 

•  The  flow  is  incompressible,  inviscid  and  axisymmetric. 

•  The  flow  is  obtained  by  superposition  of  the  mean  flow  consisting  of  a  parallel  flow  of 
velocity  U (r)  in  the  x-direction  and  an  axially  symmetric  perturbation  flow. 

•  The  mean  velocity  is  of  the  form 


IUj  if  0  <  r  <  n 

ar2  +  b  if  rx  <  r  <  r2  3  (1) 

Uoo  if  r2  <T  <  00. 

where  the  coefficients  a,  b  are  determined  so  that  U (r)  is  continuous. 

Let  Q  =  f7i  -F  q;  q  =  ui  +  wk  (i  axial,  k  radial  unit  vector).  Then 


1  The  author  is  grateful  to  Dr.  W.  W.  Bower  for  his  valuable  advice  and  support. 

2  Drazin  and  Reid:  Hydrodynamic  Stability.  Cambridge  University  Press,  (1984).  Secs,  20,21,23. 

3  This  fonn  was  chosen  for  the  axially  symmetric  case  instead  of  the  piecewise  linear  model  of 
Drazin  and  Reid,  because  "it  works". 
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Dq  — 


du  d(U  +  u)  d  (U  +  u) 

—  +  u - — - 1-  w  — 

dt  dx  dr 

dw  dui  dw 

m+uJi+wTr 


which  leads  to  the  momentum  equations: 


d 


at 


dU 


dp 


+  U—  1  U  +  W—  =  -  7T- 


dr 


dx  / 

<9  a  \  _  dp 

Ji+Udi)W~  ~  dr 


dx 


(2) 

(3) 


Continuity  equation  in  cylindrical  coordinates 


(4) 


Ansatz: 

u  =  w(r)exp[iafx  —  iut] 
p  =  p(r)exp[io:x  —  iut] 

Stream  function: 


ip(t,  x,  r )  =  ip(r)exp[iax  —  iut] 


The  factor  exp  [fax  —  iut]  will  be  assumed  and  suppressed.  The  “  ^  ”  notation 
will  be  omitted.  The  velocity  components  are  then  derivable  for  the  stream  function: 


1  d(rip) 
r  dr 


„  1  . 

=  +  ~1r, 

r 


(5) 


w  =  —  ia'tp  (6) 

Then  (6)  and  (7)  imply  the  continuity  equation. 

The  pressure  can  be  obtained  directly  from  (2).  Expressing  u,  w  in  terms  of  ip,  (2) 
can  be  rewritten  as 


„  1  , 
ip  +  -ip 
r 


where  c  =  u/a 
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Governing  equation  for  the  stream  function 

Substituting  (6)  and  (7)  into  the  momentum  equation  (3), 


(U  -  c)  o?ip  =  - 


d_ 

dr 


-  [u-^)U'  +  ^)  +  u'ip 


or 


r  -i 

r 

(U-c) 

ip" + 

r 

-{ 

\u"  +  a2(U  -  c)  -  V  +  ^(U-c)^  =  0 


Matching  conditions 

1.  Continuity  of  w(r)  : 

From  (8), 

p  d  (  ip  \  1  il> 

(U  -  cf  ~  dr\U  —  c)  rU-c 


hence 


lim 

€ — >0 


U 

and  therefore 


^r=-limr{ 

-cjr._e  e^° JTi-(  l 


p  +  lJ-Ur  =  0 


(i U-cf  rU-c 


1> 


U-c 


n+o 


=  0 


T{—  0 


and  since  U(r)  —  c  is  a  continuous  function  of  r, 


m;!2=o 


i.e.,  ip(r),  and  by  (7),  w(r)  must  be  continuous. 

2.  Continuity  of  p(r )  : 

From  Eq.  (7) 


(10') 


or,  considering  the  continuity  of  ip, 

'  d_ 
dr 

a  relation  similar  to  (6). 


V-1 

U—c 


rt+  0 

=  o 

r,-0 


Assumed  velocity  profiles:  Defined  as  in  Eq.(l), 

U (r)  is  assumed  continuous; 

Uj,  Uoo  are  known  constants; 

a,  b  are  constants  determined  from  the  continuity  of  the  velocity  profile  U  (r).  On  the 
inner  (r  =  ri)  and  outer  (r  =  r2)  boundaries  of  the  shear  layer 


ar\  +  b  =  Uj, 
ar\  +  b  — 


respectively,  from  which 

Uj  -  U0 


a  — 


r2  ~  r\ 


and  b 


r\Uj  -  r\Uo 


r2  ~  rl 


(12) 


Form  of  the  perturbation  ip(r )  in  the  three  subregions.  In  each  of  the  three 
subintervals  [0,  rx],  [ri ,  r2]  and  [r2,  oo],  it  was  assumed  that  U  (r)  is  of  the  form 

U  =  Ar2  +  B 


where  A  and  B  are  constant  in  each  of  the  intervals.  (In  fact,  A  =  0  in  the  first  and  third 
intervals).  Thus,  the  differential  equation  (8)  incredibly  simplifies  to  the  same  equation  in 
the  three  intervals: 

<  + V  +  (-«2-^)V'  =  0.  (13) 

This  is  Bessel's  differential  equation  for  the  variable  i  ar  and  parameter  value  1 .  The 
general  solution  of  (13)  is 

ip  =  CiJi(iar)  +  C^Y^iar)  =  C'iH^(iar)  +  C2  H^(iar) 


or,  alternately, 

ip  =  C\l\  (cur)  +  c2K!  (ar) 

where  I „(z)  and  Kj ,(z )  are  the  modified  Bessel  functions  defined  by 
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I„(z)  =  e-"I/2J„(e“/2z)  <14a) 

M*)  =  ■|)(ei'/sz)  =  -  Ee-i~/2H(2)(e-"/2z).  4  (14b) 

Since  i>  is  bounded  at  r  =  0  and  is  bounded  and  exponentially  decaying  as  r  -»  oo, 
and  a  is  near  real, 

f  All  (acr)  if  0  <  r  <  n 

■0(r)  =<  5Ii  (qt)  A  CKi  (or)  ifri<r<r2  (15) 

[  DKi(ar)  ifr2  <  r  <  oo 

The  functions  Ki  (z)  are  many  valued,  but  one  branch  of  each  is  real  if  z  is  real. 
These  real  branches  will  be  used  in  the  subsequent  calculations. 

Application  of  the  Matching  Conditions.  Introducing  Eqs.  (15)  into  (6)  and 
considering  (1 1)  yields 


Ali(ria)  =  5Ii(ria)  +  CK^ria)  1 

DKi(r2a)  =  Bh  (r2a)  +  CKi(r2a)  J 


(16a, b) 


Substituting  (15)  into  (30)  yields 

A(Uj-c) 


4-  +  ~  )li(<w) 

dr  r 


J  T—r  i 


=  B  [Um (r)  -  c]  +  ^li(ar)  -  2arli(ar) 


r—ri 


+  c 


[Um(r)  -  c]  (  ~  ^  )K i(or)  -  2arK1(ar) 


r=r  i 


(17a) 


4 


Encyclopedic  Dictionary  of  Mathematics  (EDM),  V.  2,  App.  A,  Table  19. IV,  Second  Ed.,  MIT 
Press,  1993;  p.  1805. 

Compare  this  with  the  only  apparently  more  precise  formuation  of  Abramowitz  Stegun  (A&S): 
Handbook  of  Mathematical  Functions  (Dover,  1965),  Eq.  9.6.4.: 


K„(z)  =  ^  -  tt  <  argz  < 

M*)  =  ~  ^  \  <  argz  <  tt^ 
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D(U oo  -  c) 


7 -  +  -Vi(ar) 

dr  r 


T=T-l 


B 


+  C 


[Um(r)-cl(£  +  1-jlI(c,r)-2arl,(ar) 


J  r=r2 


[Um(r)  -  c](  \  )Ki(ar)  -  2arKi(ar) 


J  r=7*2 


(17b) 


where  f7m(r)  =  ar 2  +  6.  The  modified  Bessel  functions  satisfy  the  recurrence  relations 

Ii  (z)  H — Ii  (^)-—  loC2) 
z 

K'i(z)+-K1(z)=  -KoW  5 

z 

Using  these  identities  we  find 

+ -^Ii  (ar)  ==  alo(ar);  + -^Kj  (ar)  =  —  aKo(ar);  (18a, b) 

hence,  one  obtains  from  (17) 


A(Uj—c) 


[i+l)h{aT) 


T~T  i 


=  B 


+  C 


[Um(r)  -  c}[  ^  ^  )li (ar)  -  2arlx(ar) 


J  r=T  i 


[C/m(r)  -  c](  7-  +  -  )  Ki  (ar )  -  2arKi(ar) 


dr  r  _ 

A(C7;— c)alo(ari) 

=  5[(Uj-c)I0(ar1)  -  2ar1Ii(ar1)] 


+  C[  -  ( Uj-c )  aKo(ari)  -  2ariKi(ari)]  J 
where  the  identity  Um{ri)  -  c  =  U— c  was  used.  Similarly, 


J  r=ri 


(19a) 


EDM,  same  location 
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-  D(Uoo  -  c)  qK0(qt2) 

=  B[(Uoo  -  c)al0(ar2)  -  2arll{ar2)}  > 

+  C[-  a(Uoo  -  c)  aK0 (ar2)  -  2ar2K1(ar2)]  ^ 

The  Dispersion  Equation 

Equations  (16a,b)  and  (19a,b)  can  be  written  in  the  concise  form 


(19  b) 


MX  =  0 


(20) 


where  X  =  [ABC  D]T  and 


li(na) 

0 

Valo(ria) 

0 


0 

(r2a) 

0 

W  aKo(o:r2) 


-  Ii(na) 

-  Ii(r2a) 

-  Valo^ri)  +  2aril1(ar1) 
W^aIo(a:r2)  -  2ar2Ii(ax2) 


-  Ki(rja) 

—  Ki(r2a) 

V«Ko(ari)  +  2ariKi(ari) 

—  W  aKo(ar2)  —  2ar2Ki(ar2) 


Here  the  abbreviations  (Uj-c)  =  V  and(C/oo  -c)  =  W  have  been  used.  The  system 
(20)  has  a  nonzero  solution  if  and  only  if  the  dispersion  equation 

A(cu,  a)  =  det  [M]  =  0  (21) 


is  satisfied.  The  matrix  M  can  be  simplified  by  elementary  row  operations.  M  is  found  to 
be  row-equivalent  to  the  matrix 


Ii  Ora) 

0 

0 

1 

rial0(ria) 

M'  = 

0 

Kj(r2a) 

1 

r2aKo(r2a) 

0 

V'alo(ria) 

0 

2oriIi  (ari) 

2ariKi(ari) 

0 

W a  Ko(qt2) 

—  2ar2I]  (ar2) 

-  2ar2Ki(ar2) 

(22) 


where  the  Wronski  determinant  identity 
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Io(2)Kx(z)  +  Ii(2)Ko(2)  —  ^ 

have  been  used.  (Abel's  theorem). 

The  dispersion  equation  detM  =  detM'  =  A(w,  a)  is  obtained  by  multiplying 
both  sides  by  rir2Q:2Io(r’itt)Ko(r2Q!),  expanding  the  determinant  and  finally,  simplifying 
by  the  same  factor. 

A(u>,  a)  =  —  rif2A(w,  a)  = 

[Uj  -  c  +  2ariIi(r1a)Ki(ria)]  [U^-c-  2ar2li(r2a)K1(r2a')]  (23) 

+  |^2arir2li(riQ:)K1(r2a:)  =0. 

where  the  coefficient  a  can  be  obtained  in  terms  of  a  from  (12),  and  c  =  u/ a. 

Note  that  det(M)  is  a  quadratic  function  of  c,  and  therefore,  one  may  express  uj 
as  a  function  of  a  explicitly.  Again,  the  dispersion  diagram  is  much  easier  to  construct  for 
the  time-unstable  case  then  for  spatial  instability.  Recommendations  about  the  numerical 
methods  available  to  solve  the  dispersion  relation  will  be  presented  at  the  end  of  this 
report. 

In  order  to  obtain  the  stream  function  ip(x,  r,  t)  corresponding  to  a  given  (tu,  a) 
pair  we  first  need  to  determine  the  coefficients  A,  B,  C,  D.  These  coefficients  are  the 
columns  of  the  4x4  matrix  M’  of  Eq.  (22)  and  satisfy  the  equation. 

M'X=0 

The  components  of  X  have  been  obtained  as  the  minors  of  M7  belonging  to  its  fourth  row 
Thus,  up  to  an  arbitrary  nonzero  constant  of  proportionality, 


(24) 


The  stream  function  ip(r)  can  be  obtained  from  Eqs.  (15a, b,c).  The  functions  u(r ), . w(r ) 
and  p(r)  are  obtained  from  (5),  (6)  and  (7),  respectively: 

6  A&S,  Eq.  9.6.15,  p.375 
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(25) 


f  Aalo(ar)  ifO  <  r  <  rq 

u(r)  =  <  Balo(ar)  +  CaK0(ar)  if n  <  r  <  r2 
[DaKo(ar)  ifr2<r<oo 

{-  iaAl\  (ar) 

—  iaBli  (ar)  —  iaCKi  (ar ) 

-  iaDKi  (ar) 


if  0  <  r  <  r\ 
if  ri  <  r  <  r2 

if  r2  <  t  <  00 


-  A(Uj  -  c)alo(ar) 


if  0  <  r  <  ri 


p(r)  =  < 


B{  -  (ar2  +  b  -  c)al0(ar)  +  2arli  (ar)] 

+  C[  -  (ar2  +  6  -  c)aKo(ar)  +  2arKi  (ar)] 


-  D(Uoo—c)aK-o(ar) 


if  ri  <  r  <  r2 


if  r2  <  r  <  00 


Proposed  numerical  approach 


(26) 


(27) 


Equations  (23),  (25),  (26)  and  (27)  define  the  solution  of  the  jet  flow  problem. 
All  these  equations  are  defined  in  terms  of  the  modified  Bessel  functions  10(2),  E  (z), 

K0  (z)  and  Ki  (z).  Therefore,  the  first  problem  is  the  evaluation  of  these  functions  for  real 
values  of  z  to  analyze  the  temporal  instability  problem,  and  for  complex  z  to  analyze 
spatial  instability.  These  functions  are  pretty  common;  in  fact,  for  real  arguments  their 
values  are  tabulated  in  Abramowitz  and  Stegun7  The  function  values  and  the  functional 
relations  can  be  also  easily  accessed  by  MATHEMATICA;  in  fact,  access  is  no  more 
difficult  than  to  the  exponential  function,  say.  The  following  table  gives  the  applicable 
notations: 


Math,  notation 

MATHEMATICA 

I  n(z) 

Bessell(n,  z) 

K  n(z) 

BesselK(n,  z) 

Alternately,  these  functions  are  not  hard  to  evaluate  in  FORTRAN.  There  is  a  subroutine 
package  available  at  Boeing,  written  by  Rick  Whitaker,  and  (to  my  understanding)  first 
developed  at  the  Argonne  National  Laboratories  in  the  1970's.  The  relevant  subroutines 
are  ZBESI  and  ZBESK.  Their  usage  is  thoroughly  documented  by  comment  lines  of  the 
source  code. 


7  A&S,  Modified  Bessel  Functions  -  Orders  0,  1  and  2  -  Table  9.8,  pp.  416-422. 
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Numerical  Solution  of  the  Dispersion  Equation 
Newton's  method 

The  application  is  straightforward:  For  the  solution  of  A  (u,  a)  =  0  for  a  fixed 
value  of  (j=ujo,  we  set  up  the  recursion 

A(an,wo) 

O’n+l  —  On  ~  ,  \  - 

AQ(an,w0) 

The  partial  derivative  in  the  denominator  is  a  simple,  but  many-term  expression 
constructed  of  the  values  of  the  modified  Bessel  functions  I  n(z),  Kn(z),  (n  =  0, 1)  and  is 
anyway  needed.  Newton's  method  will  work  well,  if  a  good  initial  guess  is  available. 

Other  iterative  algorithms,  such  as  Muller's  method,  are  also  available.  The 
latter  uses  a  quadratic  extrapolation  scheme  based  on  three  initial  points,  rather  than  the 
linear  extrapolation  of  Newton's  method  from  a  single  initial  point.  The  convergence  of 
Muller's  algorithm  is  somewhat  slower  than  that  of  Newton's.  The  advantage  of  Muller's 
method  is  that  it  does  not  require  evaluation  of  derivatives. 


Method  based  on  Gaster's  Theorem 


In  many  cases,  the  dispersion  equation  is  easier  to  solve  for  u  than  for  a.  As  a 
result,  it  is  easier  to  find  the  complex  values  of  u  for  given  real  a  than  the  complex 
solutions  for  a  if  the  real  u  is  specified.  Consequently,  the  temporal  instability  problem 
may  be  easier  to  solve  numerically  than  the  spatial  one.  This  is  indeed  the  case  for  the 
dispersion  equation  (23),  quadratic  in  terms  of  u>.  Thus,  u  —  n(a)  can  be  obtained  by 
application  of  the  quadratic  solution  formula. 

Gaster's  theorem8  allows  the  approximate  evaluation  of  the  complex  valued 
function  a  =  A(u)  if  fl(ai)  is  known  for  real  a,  provided  that  u>  and  a  are  both  small  and 
d(3a)/du  is  small  for  real  u>,  conditions  which  may  be  valid  near  neutral  stability.  Thus, 
if  the  circumstances  are  right,  the  spatial  instability  function  as  =  As(ws)  (u>s  real)  can 
be  obtained  from  the  function  uj  =  £l(aT)  describing  temporally  unstable  modes  (a  real) 
by  the  relationships 


a's 


OIt 


ft 

«s 


UJj 


(M\ 

\da'  )T 


(28) 


M.  Gaster:  A  note  on  the  relationship  between  temporally  increasing  and  spatially  increasing 
disturbances  in  hydrodynamic  stability.  JFM  Vol.  14,  pp.  222-224.  (1962) 
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where  single  primes  (')  denote  real  and  double  primes  (")  imaginary  parts.  Note  that 
a'T  =  a,  hence 

(dj\ 

\  da' )  T  da  \da )  T 

The  latter  expression  can  be  evaluated  with  reasonable  effort.  From  (24)  by  implicit 
differentiation 

fv  'N'  du 

Aa(u,  a)  +  Au(u,  a)—  =  0 


\daJT  \Aw(w,q:)Jt 


(29) 


Here,  the  denominator  is  easy  to  find,  and  the  numerator  is  a  simple  but  a  many-term 
expression  constructed  of  the  modified  Bessel  function  values  anyway  needed. 

Discussion.  The  method  outlined  gives  only  an  approximation  valid  in  the  vicinity 
of  a  =  0.  Its  result  has  to  be  refined  by  some  other  method,  such  as  Newton's.  In  other 
words,  the  result  can  serve  as  the  initial  guess  for  an  iterative  method,  and  therefore  it  is 
probably  not  our  method  of  choice. 


Tracing  a  =  A(cj)  in  the  complex  a-plane.  Analogously  to  (29),  we  may  derive  a 
differential  equation  for  a  —  A(u): 


da 

du> 


Aw{u,a) 

A«(w,o) 


(30) 


where  u  is  a  real  variable,  but  a  is  complex.  Thus,  we  have  the  system  of  real,  coupled 
differential  equations 


da'  _  ^  f  A^(a/,a//,cu)l 
du  |  Aa(a,,a",u)  J 

—  =  3m  f 

du  \  Aa(a',a",u)  J  , 


(30') 


with  the  real  unknowns  a',  a The  system  (30)  can  be  solved  numerically  with  a  Runge- 
Kutta  method  or  any  predictor-corrector  method  such  as  the  Adams-Moulton  or  Milne- 
Thompson  method.  Excellent  programs  are  available  as  canned  programs  in  packages 
such  as  the  IMSL  routines.  This  method  has  been  tested  before  by  the  author  in  the 
solution  of  the  dispersion  equation  stemming  from  another  problem.  The  method  depends 
on  the  knowledge  of  an  initial  solution  {uin.  ao}i 
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Rectangle  bisection  method. 

This  method,  developed  by  the  author  and  tested  on  dispersion  relations  from 
other  fluid  flow  problems,  is  based  on  the  connection  between  the  winding  number  and  the 
number  of  zeros  of  a  function  analytic  in  the  closure  of  an  open  domain  V  with  a  smooth 
boundary  in  the  complex  plane.  If  f(z)  is  analytic  in  the  closure  of  'D  and  the  winding 
number  of  the  map  of  the  boundary  of  V  by  /  (z)  is  u,  then  /  ( z )  has  v  zeros  in  T> . 
(Argument  principle) 

Here  only  a  sketch  of  the  method  can  be  given.  The  argument  principle  is  applied 
to  a  rectangular  domain  H  of  the  a-plane  which  does  not  contain  any  singularity  or  branch 
cut  of  the  function  A  (c o0,  a)!  Suppose  that  by  the  argument  principle  it  is  determined 
that  H  contains  a  positive  number  of  zeros  of  A  {uJ0,a).  Then,  H  is  divided  by  a 
horizontal  (or  vertical)  center  line  into  the  rectangles  Hi  and  H2.  Apply  the  argument 
principle  to  Hi  and  H2.  If  the  function  values  of  /  on  the  boundary  of  H  are  saved,  we 
need  to  determine  only  the  function  values  on  the  center  line  segment  used  in  the  sepration 
o{Hx  and  H2.  The  distribution  of  the  sampling  points  for  the  evaluation  of  /  is  optimized 
by  an  adaptive  method.  Roughly,  the  step  size  of  function  evaluations  is  halved  and 
doubled  if  the  argument  variation  gets  too  large  and  too  small,  respectively.  If  the 
rectangle  Hi  contains  a  positive  number  of  zeros  we  save  the  coordinates  of  Hi  and  the 
number  of  zeros  it  contains,  otherwise  we  discard  the  information.  The  rectangle  H2  is 
treated  the  same  way.  Now  the  algorithm  is  repeated  on  each  of  the  remaining  rectangles 
involving  bisection  of  the  rectangles,  finding  numbers  of  zeros  contained  in  the  new 
rectangles,  and  discarding  any  empties.  The  direction  of  the  bisection  is,  in  general, 
alternating  between  horizontal  and  vertical.  We  continue  the  iteration  until  the  rectangle 
diagonals  are  smaller  than  a  specified  error  tolerance.  (See  Fig.  1). 

The  method  has  been  coded  by  the  author  and  applied  to  several  unrelated 
problems  The  method  is  always  convergent  and  convergence  is  surprisingly  fast  (though, 
of  course  not  comparable  to  the  speed  of  convergence  of  Newton's  algorithm,  if  it 
converges).  The  speed  of  convergence  depends  on  the  time  required  for  a  single  function 
evaluation.  Roughly,  each  new  “generation”  of  rectangles  requires  4  to  10  function 
evaluations  and  the  number  of  generations  can  be  obtained  from  the  error  formula 


where  do  is  the  initial  rectangle  size  and  n  is  the  number  of  generations. 

Discussion.  It  was  found  in  the  planar  case  that  if  a  is  plotted  in  the  complex  en¬ 
plane  and  the  real  cu  is  a  parameter,  the  resulting  dispersion  diagram  has  infinitely  many 
branches  starting  from  points  a°k  corresponding  to  u  =  0.  It  was  found  that  a0  =  0,  but 
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fRc[ak}~*  -  oo9and3m[afc]->ooas  k-*oo.  (See  Fig.  2)  Clearly,  Gaster's  theorem 
cannot  be  used  to  find  initial  points  on  the  higher  branches.  On  each  branch,  a  starting 
point  is  needed.  This  can  be  obtained  by  a  “lucky  guess”,  or  a  few  steps  of  the  rectangle 
bisection  method.  Whichever  method  we  used  to  find  an  initial  guess,  it  must  be  followed 
by  refinement  by  Newton's  method,  or  some  other  similar  iterative  algorithm. 


The  calculation  was  based  on  the  dispersion  relation 

u2  =  (1  -  2a)2  -  e~4a 

of  Drazin  and  Reid.  They  base  this  on  the  assumption  that  [/^  =  +  1  and  £/_ «,  =  -  1,  hence, 
Uoo  >  U- oo.  In  our  present  model  the  opposite  is  the  case.  This  accounts  for  the  sign  difference 
of  £He  [a] 
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Evolution  strategies  for  parameter 
optimization  in  jet  flow  control 

By  P.  Koumoutsakos,  J.  Freund1  AND  D.  Parekh2 


We  present  results  from  the  application  of  evolution  strategies  for  parameter  opti¬ 
mization  in  direct  numerical  simulations  and  vortex  models  of  controlled  jet  flows. 
It  is  shown  that  evolution  strategies  are  a  portable,  highly  parallel  method  that  can 
complement  our  physical  intuition  in  the  parameter  optimization  of  such  flows. 


1.  Introduction 

For  centuries  engineers  have  taken  inspiration  from  nature  in  designing  efficient 
aerodynamic  configurations.  It  is  no  coincidence  that  the  shape  of  an  aircraft  s 
wing  resembles  a  bird’s.  We  wish  to  approach  the  problem  of  flow  control,  not 
from  the  perspective  of  imitating  existing  natural  forms,  but  from  the  perspective 
of  developing  efficient  control  algorithms,  by  employing  techniques  inspired  by  bi¬ 
ological  processes.  These  techniques,  which  we  will  refer  to  as  “machine  learning 
algorithms” ,  are  gaining  significance  in  the  areas  of  modeling  and  optimization  for 
fluid  dynamics  problems  as  a  technology  that  could  help  reduce  cost  and  time  to 
market  of  new  designs. 

1.1  Evolution  strategies 

Some  of  the  seminal  work  in  this  field  (Rechenberg  1971,  Schwefel  1974,  Hoffmeis- 
ter  1991)  actually  was  aimed  at  improving  aerodynamic  shapes.  As  stated  in 
(Schwefel,  1974): 

“In  1963  two  students  at  the  Technical  University  of  Berlin  met  and  were  soon 
collaborating  on  experiments  which  used  the  wind  tunnel  of  the  Institute  of  Flow 
Engineering.  During  the  search  for  the  optimal  shape  of  bodies  in  a  flow,  which 
was  then  a  matter  of  laborious  intuitive  experimentation,  the  idea  was  conceived  of 
proceeding  strategically.  However,  attempts  with  the  coordinate  and  simple  gradi¬ 
ent  strategies  were  unsuccessful.  Then  one  of  the  students,  Ingo  Rechenberg,  now 
professor  of  Bionics  and  Evolutionary  Engineering,  hit  upon . the  idea  of  trying  ran¬ 
dom  changes  in  the  parameters  defining  the  shape,  following  the  example  of  natural 
mutations.  The  evolution  strategy  was  bom.  ”  (The  second  student  was  Hans  Paul 
Schwefel) . 

Since  this  pioneering  work,  stochastic  optimization  techniques  have  gained  recog¬ 
nition  and  popularity  in  several  fields  of  engineering,  but  this  has  not  been  the  case 
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in  the  field  of  fluid  dynamics  in  the  last  three  decades.  Recent  work  by  Rechen- 
berg  (1994)  focuses  on  the  shape  optimization  approach  with  the  construction  and 
experimental  testing  of  shapes  that  have  been  produced  via  evolutionary  strate¬ 
gies  using  computer  simulations.  Evolution  strategies  have  also  been  implemented 
in  order  to  optimize  the  motions  of  an  artificial  tuna  (M.  Triantafvllou,  private 
communication) . 

Here,  we  report  preliminary  results  from  the  application  of  evolution  strategies 
in  the  optimization  of  actuator  parameters  in  active  jet  flow  control  and  in  the 
optimization  of  bifurcating  and  blooming  jets. 

1.2  Jet  flow  control 

It  is  desirable  in  many  circumstances  to  enhance  mixing  in  the  exhaust  from 
aircraft  engines.  Applications  include  lift  enhancement,  signature  reduction,  and 
temperature  reduction  on  blown  flaps.  This  work  focuses  on  the  latter  case.  The 
blown  flap  on  a  C-17  (Fig.  1)  is  currently  made  out  of  titanium  to  avoid  melting. 
If  mixing  can  be  significantly  enhanced  so  that  the  plume  temperature  is  reduced, 
the  flap  could  be  constructed  from  aluminum,  a  much  less  heavier  and  expensive 
alternative. 


Figure  1.  Blown  flap  as  on  a  C-17. 

Recently,  actuators  have  been  developed  and  tested  on  a  full-scale  engine  which 
have  the  control  authority  to  accomplish  this  objective.  The  goal  of  this  work  is  to 
optimize  their  parameters  to  maximize  their  effectiveness.  This  is  being  undertaken 
as  a  joint  experimental,  numerical,  and  control  theory  effort.  The  discussion  here  is 
limited  to  the  simulations  and  the  application  of  evolution  strategies  to  the  problem. 

1.3  Optimization  of  bifurcating  and  blooming  jets 

The  proper  combination  of  axial  and  helical  excitation  at  different  frequencies 
generates  the  unique  class  of  flows  known  as  bifurcating  and  blooming  jets  (Lee  and 
Reynolds  1985,  Parekh  et  al.  1987).  The  axial  forcing  causes  the  shear  layer  to  roll 
up  into  distinct  vortex  rings  at  the  forcing  frequency.  The  helical  excitation  perturbs 
the  rings  radially,  producing  a  small  eccentricity  in  the  ring  alignment.  This  initial 
eccentricity  is  amplified  by  the  mutual  ring  interactions  leading  to  dramatic  changes 
in  jet  evolution.  When  the  axial  frequency  is  exactly  twice  that  of  the  helical 
excitation,  the  jet  bifurcates  into  two  distinct  jets,  with  successive  rings  moving 
alternately  on  one  of  two  separate  trajectories.  This  Y-shaped  jet  spreads  at  angles 
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over  80  degrees,  depending  on  forcing  frequency  and  amplitude.  The  relative  phase. 
<p.  between  the  axial  and  helical  forcing  signals  determines  the  plane  in  which  the 
jet  bifurcates.  When  the  ratio,  0,  of  axial  to  helical  excitation  frequency  is  non¬ 
integer,  the  vortex  rings  scatter  along  a  conical  trajectory.  When  viewed  from 
downstream,  the  vortex  ring  pattern  often  resembles  a  bouquet  of  flowers,  hence 
the  name  “blooming  jet.” 

In  applying  the  evolution  strategies  to  this  class  of  flows,  we  are  exploring  whether 
the  phenomena  discovered  experimentally  could  also  be  obtained  in  our  simulations 
via  an  “evolutionary  process”  and  whether  new  phenomena  could  be  found.  Here 
a  vortex  model  describes  the  jet  dynamics.  The  optimization  algorithm  is  tuned  to 
maximize  jet  spreading  by  varying  the  excitation  parameters. 

2.  Evolution  strategies  for  optimization 

We  discuss  first  the  formulation  of  evolution  strategies  for  the  optimization  of 
N-dimensional  functions: 


F(x)  =  F(x1,x2,...,xm) 

We  define  a  vector  in  the  parameter  space  as  an  individual.  The  whole  discrete 
parameter  space  can  then  be  considered  as  a  population  of  individuals.  Evolution 
strategies  try  to  identify  the  best  individual  from  this  population  based  on  the  fitness 
value,  prescribed  by  the  function  F.  The  optimization  proceeds  by  following  to  a 
certain  extent  models  of  biological  evolution. 

2.1  Two  membered  evolution  strategies 

The  simplest  (and  earliest)  form  of  evolution  strategies  is  based  on  populations 
that  consist  of  two  competing  individuals  ( “a  two-membered  strategy”).  The  evo¬ 
lution  process  consists  of  the  two  operations  that  Darwin  (1859)  considered  as  the 
most  important  in  natural  evolution:  mutation  and  selection.  Each  individual  (i.e. 
vector  in  the  parameter  space)  is  represented  using  a  pair  of  floating  point  valued 
vectors: 

u  =  u(x,  cr) 

where  cr  is  an  M-dimensional  vector  of  standard  deviations. 

Following  Rechenberg  (1971)  and  using  terminology  from  biology,  the  optimiza¬ 
tion  algorithm  may  be  described  as  follows: 

a  -  Initialization:  A  parent  genotype  consisting  of  M- genes  is  specified  initially  (x°). 
At  each  generation  an  individual  u £  =  (x£ ,  er™ )  is  identified. 

b  -  Mutation:  The  parent  of  generation-n  produces  a  descendant,  whose  genotype 
differs  slightly  from  that  of  the  parent.  The  operation  of  mutation  is  then  realized 
by  modifying  x  according  to: 


xc  =  xp  +  A/’(0,<Xp)  2.1 

where  J\f(0,  cr)  denotes  an  M-dimensional  vector  of  random  Gaussian  numbers 
with  zero  mean  and  standard  deviations  cr. 
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c  -  Selection:  Due  to  their  different  genotypes  the  two  individuals  of  the  population 
can  have  a  different  fitness  for  survival.  This  fitness  is  evaluated  by  the  function 
/.  Only  the  fittest  of  the  two  individuals  is  allowed  to  produce  descendants  at 
the  following  generation.  Hence  to  minimize  F  we  write: 

Xn+1  =  f  x£,  if  F(x»)  <  F(x-);  2  2 

p  \xc>  otherwise. 

Note  that  in  this  two-membered  algorithm  the  vector  cr  of  standard  deviations 
remains  unchanged  throughout  the  evolutionary  process. 

For  regular  optimization  problems  (see  Michalewicz,  1996  for  a  definition)  it  is 
possible  to  prove  the  convergence  of  the  method  to  a  global  minimum.  However, 
this  theorem  does  not  provide  a  convergence  rate  of  the  method. 

In  this  work  we  have  implemented  the  1/5  success  rule  proposed  by  Rechen- 
berg  (1971).  According  to  this  rule:  During  the  optimum  search  the  frequency 
of  successful  mutations  is  checked  periodically  by  counting  the  ratio  of  the  number 
of  successes  to  the  total  number  of  trials.  The  variance  is  increased  if  this  ratio  is 
greater  than  1/5  and  it  is  decreased  if  it  is  less  than  1/5.  The  period  over  which  this 
performance  is  being  checked  depends  on  the  number  of  parameters  that  are  being 
optimized.  We  refer  to  Schwefel  (1995)  for  further  details  on  the  implementation  of 
the  two-membered  evolution  strategies. 

2.2  Multi-membered  evolution  strategies 

One  of  the  drawbacks  of  the  1/5  rule  for  the  two-membered  strategy  is  that  it 
may  lead  to  premature  convergence,  as  the  step  lengths  can  be  reduced  to  zero,  thus 
not  improving  the  progress  towards  a  global  optimum.  There  are  several  possible 
remedies  to  this  drawback.  Of  particular  interest  are  those  that  can  be  constructed 
by  further  developing  the  model  of  evolution  to  resemble  natural  processes.  In  that 
context,  a  higher  level  of  imitation  of  an  evolutionary  process  can  be  achieved  by 
increasing  the  number  of  members  in  a  population.  Such  multi-membered  strategies 
are  usually  formulated  in  terms  of  p- parents  and  A-descendants.  The  most  common 
strategies  are  then  described  as  (/r,A)  and  (p  +  A).  In  the  (p,  A)  case  at  each 
generation  the  p-fittest  individuals  are  selected  only  among  the  A  children  of  the 
generation,  whereas  in  the  (p+ A)  case  the  parents  are  also  included  in  the  evaluation 
process.  Schwefel  (1995)  presents  an  extensive  comparison  of  multi-membered  and 
two-membered  evolution  strategies  for  a  series  of  optimization  problems. 

2.3  Handling  of  constraints 

One  of  the  advantages  of  evolution  strategies  is  the  ease  and  simplicity  by  which 
they  can  handle  problem  constraints.  Such  constraints  are  usually  formulated  as 
inequalities.  For  example  in  the  case  of  q  constraints  of  the  parameters  x  we  require 
that: 

Cj(x)  >  0  for  all  j  =  1,  ...,q 

Descendants  of  a  certain  parent  that  do  not  satisfy  the  constraints  are  accounted  as 
results  of  unsuccessful  mutations.  Occasionally  the  boundaries  of  the  constrained 
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FIGURE  2.  Schematic  shows  the  nozzle  (a)  and  actuators  (b). 

regions  are  smoothed  out  in  order  to  facilitate  the  convergence  of  the  method  in 
highly  constrained  problems. 

3.  Jet  flow  control 

The  compressible  flow  equations  were  solved  with  direct  numerical  simulation 
using  a  combination  of  sixth  order  compact  finite  differences,  spectral  methods, 
and  fourth  order  Runge  Kutta  time  advancement.  Further  details  of  the  numerical 
algorithm  and  techniques  for  including  actuators  into  the  calculations  were  recently 
reported  by  Freund  &  Moin  (1998).  Naturally,  in  a  direct  numerical  simulation  we 
are  restricted  to  highly  simplified  geometries  (Fig.  2);  nevertheless,  the  actuators 
were  able  to  reproduce  the  effects  observed  in  experiments  by  Parekh  et  al.  (1996). 
Figure  3  shows  a  visualization  of  a  jet  forced  into  a  flapping  mode  and  an  unforced 
jet.  Clearly,  the  mixing  is  enhanced  downstream. 

For  this  preliminary  study,  only  three  types  actuation  parameters  were  varied: 
the  amplitude,  frequency,  and  phase.  The  actuation,  was  a  simple  waveform  sum  of 
harmonic  waveforms: 

Vr  =  XJ Ai  +  sin  (~JTt  +  ssn(cos^))  »  3,1 

where  vr  is  the  radial  velocity  at  the  actuator  exit  and  A{  are  the  amplitudes,  St{  are 
the  Strouhal  numbers,  and  fa  are  the  phases  of  the  different  modes.  The  sgn(cos(#)) 
causes  each  waveform  to  excite  a  flapping  mode  in  the  jet.  Note  that  the  phases, 
fa,  are  the  relative  phases  of  the  different  modes  setting  the  two  actuators  always 
at  180°  out  of  phase.  The  flow  rate  out  of  either  actuator  was  constrained  to  be  less 
than  U/2  where  U  is  the  jet  velocity.  This  was  accomplished  by  simply  “clipping” 
the  velocities  to  be  below  this  level. 

The  only  constraint  on  A{  was  that  they  be  non-negative.  Strouhal  numbers  were 
restricted  to  be  0  <  St  <  0.8  and  the  phases  were  constrained  to  be  fa  €  [0,  2tt\. 

A  very  low  Reynolds  number  (Re  =  500)  jet  at  Mach  0.9  was  simulated  in  this 
preliminary  effort  to  minimize  the  computational  expense.  The  computational  mesh 
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FIGURE  3A.  Unforced  turbulent  jet.  Visualization  of  vorticity  magnitude. 


was  112  x  42  x  16  in  the  streamwise,  radial,  and  axial  direction  respectively  and  the 
computational  domain  extended  to  16  radii  downstream  and  5  radii  in  the  radial 
direction.  A  stretched-mesh  boundary  zone  was  positioned  outside  of  the  region  to 
cleanly  absorb  fluctuations  convecting  out  of  the  domain.  In  each  iteration  of  the 
evolution  strategy,  the  jet  was  simulated  starting  from  an  unforced  case  for  several 
periods  of  forcing  after  the  passing  of  initial  transients.  Because  the  flow  becomes 
quasi-periodic,  this  was  sufficient  to  provide  a  measure  of  the  long-time  actuator 
effectiveness.  Each  iteration  required  approximately  10  minutes  and  in  total  200 
iterations  were  made  (the  best  case  was  found  after  approximately  150  iterations). 

Three  wave  forms  ( N  =  3)  where  used  and  the  initial  control  parameters  were 
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5  '  ’  To 


FIGURE  4 A. Jet  mixture  fraction 
for  the  first  guess  parameters. 


FIGURE  4B.  Jet  mixture  fraction 
with  the  best  case  parameters  af¬ 
ter  200  iterations. 
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The  goal  set  for  the  evolution  strategy  was  to  maximize 

poo  p2n  pSrQ 

Q=  /  /  v^rdrdddx. 
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This  metric  Q,  was  increased  by  over  a  factor  of  10  from  the  initial  case  by  the  best 
case  parameters: 
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It  is  interesting  to  note  that  the  evolution  strategy  “chose”  to  reduce  the  amplitude 
of  two  of  the  wave  modes  to  a  very  low  level.  Effectively,  it  found  the  same  ad 
hoc  scheme  that  was  shown  to  be  successful  by  Parekh  et  al.  (1996)  and  Freund  & 
Moin  (1998).  A  forced  and  unforced  case  are  visualized  in  Fig.  4.  The  best  case 
clearly  shows  a  high  amplitude  flapping  mode  which  would  greatly  enhance  mixing 
downstream. 

4.  Vortex  model  of  bifurcating  and  blooming  jets 

In  this  work  we  model  a  circular  jet  by  the  combination  of  discrete  vortex  filaments 
and  a  semi-infinite  cylindrical  sheet  of  vorticity.  The  cylindrical  sheet  models  the 
nozzle  source  flow  whereas  the  ring  filaments  model  the  vortex  rings  generated  by 
the  axial  excitation  of  the  shear  layer. 

The  semi-infinite  sheet  of  vorticity  extends  from  —  oo  to  the  origin.  Its  axis  defines 
the  jet  centerline,  and  the  end  of  the  sheet  is  identified  with  the  jet  exit.  The 
helical  excitation  used  in  the  experiments  of  Lee  and  Reynolds  (1985)  is  modeled 
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by  rotating  the  axis  of  the  vortex  cylinder  about  the  nominal  jet  centerline.  The 
displacement,  A h,  of  the  jet  centerline  from  the  nominal  centerline  corresponds  to 
the  amplitude  of  excitation,  and  Ah  =  Ah/R.  The  rotation  frequency  is  given  by: 

A  =  4.1 

Kf 


where  the  orbital  frequency  is  defined  as: 


f*=St 4-2 

The  frequency  fa  is  the  rate  at  which  filaments  are  generated  at  the  origin. 

The  interaction  of  the  vortex  sheet  with  the  filaments  is  assumed  to  be  such  that 
the  sheet  influences  the  motion  of  the  filaments  but  the  filaments  do  not  influence 
the  sheet.  The  velocities  induced  by  each  filament  and  by  the  jet  function  are 
superimposed  to  determine  the  trajectory  of  each  filament.  The  Strouhal  number 
sets  the  time  between  creation  of  new  ring  filaments  at  the  origin. 

The  circulation  of  each  filament  is  identical  and  is  determined  from  circulation 
conservation  constraints.  Assuming  the  thickness  of  the  cylindrical  sheet  to  be 
much  smaller  than  its  radius,  the  vorticity  flux  (per  unit  of  circumference)  within 
the  sheet  through  any  plane  perpendicular  to  the  jet’s  axis  is  given  by  U2/2.  By  the 
assumption  of  a  perfect  fluid,  the  vorticity  convected  from  the  cylindrical  sheet  must 
equal  the  vorticity  convected  by  the  discrete  filaments.  This  conservation  relation 
can  be  expressed  in  terms  of  T  and  7  as 


where  T  is  the  circulation  of  each  ring  filament,  7  is  the  circulation  per  unit  length  of 
the  cylindrical  vortex  sheet,  and  At  is  the  time  between  generation  of  ring  filaments. 
By  Eqs.  4.8  and  4.9,  one  obtains 


7  =  St„|.  4-4 

Further  details  concerning  the  applicability  of  this  model  and  its  numerical  im¬ 
plementation  axe  reported  in  Parekh  et  al.  (1988). 

4.1  Parameter  optimization  using  evolution  strategies 

The  primary  parameters  that  govern  the  jet  evolution  Sta,  (3  (frequency  ratio  of 
axial  and  orbital  excitation),  Aa,  Ah,  and  (f>.  The  effect  of  the  axial  excitation,  Aa. 
is  approximated  by  generating  distinct  vortex  rings  at  the  axial  forcing  frequency. 
The  sensitivity  to  axial  forcing  amplitude  is  not  modeled.  In  these  simulations  the 
other  four  parameters  are  allowed  to  vary  over  the  following  ranges:  0  <  Aa  <  1, 
0.1  <  Sta  <  1,  0.2  <  (3  <  5,  0  <  (j)  <  2n.  Different  flow  patterns  can  be  observed 
with  variations  in  (3  for  fixed  values  of  the  other  parameters.  The  simulation  is  able 
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FIGURE  4. 1A.  Blooming  jets:  Ex¬ 
perimental  Results  of  Lee  and  Reynolds  FIGURE  4.  IB.  Blooming  Jets:  Sim- 

(1985) :  f3  =  1.7,  Ah  =  0.04,  Sta  =  ulations  /3  =  1.7,  Ah  =  0.05,  Sta  = 

0.46,  Re  =  4300.  °-55- 

to  represent  qualitatively  the  full  range  of  jet  phenomena  observed  in  experiments, 
including  bifurcating  and  blooming  jets  (Fig.  4.1). 

For  the  optimization,  several  metrics  for  jet  spreading  angle  were  considered, 
including  the  average  radial  displacement  of  the  vortex  elements,  jet  spreading  angle, 
and  ring  trajectory  angles.  We  also  considered  amplitude  normalized  formulations 
of  these  metrics  to  account  for  the  cost  of  excitation.  The  metrics  were  evaluated 
over  a  broad  range  of  test  cases  to  check  if  they  would  be  robust  enough  to  provide 
the  proper  relative  rating  over  the  parameter  space  considered.  Some  metrics  are 
artificially  biased  by  the  initial  displacement  of  the  rings  or  by  normalization  with 
very  small  excitation  amplitudes.  One  metric  that  is  both  simple  and  effective  for 
this  simulation  is  the  average  angle  of  the  nominal  ring  trajectories.  For  each  case, 
this  metric  is  evaluated  after  the  same  number  of  periods  (typically,  eleven)  of  axial 
excitation.  The  nominal  ring  trajectory  angle,  9,  is  defined  as  the  angle  between 
the  jet  centerline  and  the  line  that  connects  the  center  of  the  jet  exit  to  the  centroid 
of  the  vortex  ring  nodes. 

Starting  with  an  initial  guess  for  each  of  these  parameters  and  constraints  on 
the  range  of  values  allowable  for  each  parameter,  the  genetic  algorithm  searches 
to  optimize  jet  spreading.  The  scope  of  this  work  did  not  allow  for  an  exhaustive 
investigation  of  the  parameter  space  and  convergence  characteristics,  but  even  these 
preliminary  simulations  yielded  promising  results.  With  all  four  parameters  varied 
simultaneously,  the  genetic  algorithm  selects  a  blooming  jet  similar  to  what  has 
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FIGURE  4.2A.  Hybrid  bifurcating 
jet  with  Sta  =  0.28,  Ah  =  0.63,  / 3  = 
2, and$  =  0  (side  view). 


FIGURE  4.2B.  End  view.  Each 
ring’s  32  nodes  are  plotted  as  solid 
circle. 


been  observed  in  experiments. 

The  most  striking  result  was  found  when  we  constrained  (3  =  2  and  kept  (j)  fixed. 
Initially  we  expected  the  algorithm  to  select  a  bifurcating  jet  similar  to  Fig.  4.1  A 
with  values  of  StanndAh  that  maximize  the  spreading  angle.  Instead,  a  unique  jet 
flow  (Fig.  4.2)  was  found  that  had  never  been  observed  in  previous  experiments  or 
calculations.  This  jet  flow  initially  resembles  a  bifurcating  jet.  Several  diameters 
downstream,  however,  the  two  branches  of  the  jet  exhibit  a  secondary  bifurcation 
in  which  the  rings  change  direction  along  a  path  with  an  azimuthal  angle  about 
7r /4  different  from  their  original  trajectory.  This  results  in  a  wide  spreading  angle 
as  seen  in  Fig.  4.2B. 

The  simulation  often  has  difficulty  providing  valid  solutions  for  Sta  >0.4  since  the 
initial  ring  filaments  get  tangled  together  and  quickly  degrade  to  an  unrealistic  state. 
This  constraints  were  implemented  in  the  evolution  strategy  by  simply  considering 
these  cases  as  unsuccessful  tries  for  the  optimization  algorithm. 


5.  Summary  and  conclusions 

These  preliminary  results  from  the  application  of  evolution  strategies  to  the  prob¬ 
lem  of  flow  control  suggest  that  stochastic  optimization  can  be  a  valuable  tool  that 
can  complement  physical  understanding  and  deterministic  optimization  techniques. 
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As  a  closing  remark,  we  quote  from  Schwefel: 

Since  according  to  the  'No- Free- Lunch ”  (NFL)  theorem  (Wolpert  and  Macready. 
1996)  there  cannot  exist  any  algorithm  for  solving  all  optimization  problems  that 
is  on  average  superior  to  any  competitor,  the  question  of  whether  evolutionary  al¬ 
gorithms  are  inferior /superior  to  any  alternative  approach  is  senseless.  The  NFL 
theorem  can  be  corroborated  in  the  case  of  EA  versus  many  classical  optimization 
methods  insofar  as  the  latter  are  more  efficient  in  solving  linear,  quadratic,  strongly 
convex,  unimodal,  separable,  and  many  other  problems.  On  the  other  hand,  EA 's 
do  not  give  up  so  early  when  discontinuous,  nondifferentiable,  multimodal,  noisy, 
and  otherwise  unconventional  response  surfaces  are  involved.  Their  robustness  thus 
extends  to  a  broader  field  of  applications,  of  course  with  a  corresponding  loss  of 
efficiency  when  applied  to  the  classes  of  simple  problems  classical  procedures  have 
been  specifically  devised  for.  ” 

Hence,  in  the  realm  of  flow  control,  the  key  issue  is  the  identification  of  a  suitable 
optimization  method  for  the  specific  problem  in  hand.  The  portability,  ease  of 
parallelization,  and  the  results  reported  herein  and  in  (Muller  et  al.  1999),  suggest 
that  EA’s  present  a  powerful  technique  for  parameter  optimization  in  problems  of 
flow  control. 
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Appendix  F 

Key  Results  from  Small-Scale  Engine  Experiment 

This  appendix  provides  summary  plots  representative  of  the  key  results  from  the 
experimental  demonstration  on  the  small-scale  engine.  Figure  F-l  presents  the  centerline 
temperature  decay  of  the  baseline  flow  at  two  different  Mach  numbers.  The  axial 
distance  is  normalized  by  the  nozzle  exit  diameter.  Figures  F-2  through  F-4  provide  the 
corresponding  Strouhal  number  dependence  of  centerline  temperature  at  a  fixed  location 
of  X/D  =  4.  From  a  feedback  control  standpoint,  one  would  like  to  be  able  to  average 
over  a  short  time  when  measuring  the  temperature  sensor  data  to  shorten  the  feedback 
loop  time.  This,  however,  results  in  significant  uncertainty  and  variability  in  the 
information  being  fed  back  to  the  control  loop  (see  Figures  F-3  and  F-4).  Consequently, 
a  gradient  search  technique  would  be  impractical  to  implement  because  of  the  difficulty 
of  computing  meaningful  derivatives  from  the  measured  discrete  data.  In  contrast,  the 
genetic  algorithms  are  able  to  handle  this  uncertainty  with  little  difficulty. 

Optimization  experiments  are  presented  in  Figures  F-5  through  F-7  for  a  pair  of  Mach 
numbers  and  varying  initial  guesses.  Comparison  of  different  experiments  initiated  with 
an  identical  guess  (Figure  F-5)  illustrates  the  variability  of  the  convergence  due  to  the 
inherent  randomness  in  the  process.  Figure  F-6  compares  optimization  when  the  initial 
guess  is  close  to  the  global  or  local  minimum.  The  simple  approach  implemented  here 
can  get  trapped  in  a  local  minimum  for  a  period  of  time;  with  more  iterations,  the 
algorithm  would  eventually  find  the  global  minimum.  Figure  F-7  provides  a 
representative  optimization  experiment  at  a  higher  Mach  number  of  0.32. 

Finally,  we  explore  the  robustness  of  the  technique  to  changes  in  the  system  that  are 
not  known  by  the  algorithm.  This  is  simulated  by  varying  the  engine  speed  several  times 
during  a  single  optimization  experiment.  No  information  about  the  engine  speed  nor 
even  the  fact  that  a  change  occurred  is  monitored  by  the  controller.  Without  any 
knowledge  of  the  plant,  the  nature  of  the  algorithm  enables  it  inherently  to  adapt.  While 
much  can  be  done  to  improve  the  efficiency  of  this  adaptation,  Figure  F-8  demonstrates 
the  ability  of  even  this  simple  algorithm  to  adapt  to  major  changes  in  the  plant. 
Suitability  of  this  adaptation  for  real-time  control  in  an  operational  mode  (in  contrast  to 
real-time  optimization  in  a  test  mode)  would  depend  on  the  acceptability  of  running  the 
actuation  at  random  points  in  the  parameter  space  during  operations. 
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Average  T*  Values  for  Characterization  of  Strouhal  Dependance  (Ma  =  0.15,  5%  MFR) 


Figure  F-3.  Data  scatter  for  short-time  (0.5-ls)  average  measurements  of  T*  variation  with 
temperature  (upper  plot)  and  the  corresponding  average  profile  and  polynomial  Fit  (lower 
plot)  for  M=  0.15  with  5%  mass  injection,  measured  at  X/D  =  4. 
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Genetic  Algorithm  Convergence  (Ma  =  0.15,  5%  MFR,  10  Generations/cycle) 
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Genetic  Algorithm  Convergence  (Ma  =  0.15,  5%  MFR,  10  Generations/cycle) 
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Figure  F-5.  Comparison  of  two  optimization  experiments  run  at  the  same  conditions  with 
the  identical  initial  guess  for  Strouhal  number. 


Genetic  Algorithm  Convergence  when  Initial  Guess  is  close  to  Global  Minimum 
(Ma  =  0.15,  5%  MFR,  10  Generations/cycle) 
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Genetic  Algorithm  Convergence  when  Initial  Guess  is  close  to  Local  Minimum 
(Ma  =  0.15,  %MFR  »  5%,  10  Generations/cycle) 
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Figure  F-6.  Comparison  of  optimization  experiments  run  at  the  same  conditions  but  with 
differing  initial  guesses:  close  to  global  minimum  (upper  plot)  and  close  to  local  minimum 

(lower  plot). 


113 


Genetic  Algorithm  Convergence  {Ma  =  0.32,  2.5%  MFR,  10  Generations/cycle) 
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Figure  F-7.  Example  of  optimization  at  M=0.32. 
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Genetic  Algorithm  Convergence  for  Changing  Engine  Speeds 
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Genetic  Algorithm  Convergence  for  Changing  Engine  Speeds 


Figure  F-8.  Optimization  experiment  in  the  presence  of  changing  system  conditions. 
Engine  speed  is  varied  at  various  points  during  the  optimization  process  without  any 
feedback  to  the  algorithm  that  the  engine  speed  has  been  changed. 
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